Octuple-precision

.github/workflows/codecov.yml

@@ -1,5 +1,10 @@
 steps:
 - uses: actions/checkout@master
+- name: Generate Coverage Data
+  run: |
+    JULIA_COVERAGE=true julia runtests.jl
 - uses: codecov/codecov-action@v3
   with:
     token: ${{ secrets.CODECOV_TOKEN }}
+    commit: HEAD~5..HEAD # Analyze the last 5 commits
+

.github/workflows/documentation.yml

@@ -14,7 +14,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.4'
+          version: '1.9.2' # Specify the required Julia version here
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

Project.toml

@@ -1,7 +1,7 @@
 name = "sparseQFCA"
 uuid = "077cf62c-ad2d-5edd-99ec-638372f8b004"
 authors = ["Mojtaba Tefagh <mtefagh@sharif.edu>", "Iman Ghadimi <iman.diligent@gmail.com>"]
-version = "1.9.0"
+version = "2.0.0"
 
 [deps]
 COBREXA = "babc4406-5200-4a30-9033-bf5ae714c842"
@@ -14,10 +14,15 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
+CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 
 [compat]
-COBREXA = "1.0, 1.2"
-Colors = "0.12"
-GLPK = "0.14, 1.0"
-JuMP = "0.21, 1.0"
-julia = "1.8, 1.9"
+COBREXA = "1.5.1"
+MathOptInterface = "1.28.1"
+GLPK = "1.1.3"
+Clarabel = "0.7.1"
+Colors = "0.12.10"
+JuMP = "1.21.1"
+julia = "1.9.2"

src/ConsistencyChecking/SwiftCC.jl

@@ -10,7 +10,9 @@ module SwiftCC
 
 export Model_CC, model_CC_Constructor, swiftCC
 
-using GLPK, JuMP, COBREXA, LinearAlgebra, SparseArrays, Distributed
+using GLPK, JuMP, COBREXA, LinearAlgebra, SparseArrays, Distributed, Clarabel
+
+import CDDLib
 
 include("../Pre_Processing/Pre_processing.jl")
 
@@ -66,7 +68,7 @@ and assigns values to its fields based on the other arguments passed in.
 """
 
 function model_CC_Constructor(ModelObject_CC::Model_CC, S::Union{SparseMatrixCSC{Float64,Int64}, AbstractMatrix}, Metabolites::Array{String,1},
-                              Reactions::Array{String,1}, Genes::Array{String,1},m::Int, n::Int, lb::Array{Float64,1}, ub::Array{Float64,1})
+                              Reactions::Array{String,1}, Genes::Array{String,1}, m::Int, n::Int, lb::Array{Float64,1}, ub::Array{Float64,1})
      ModelObject_CC.S = S
      ModelObject_CC.Metabolites = Metabolites
      ModelObject_CC.Reactions = Reactions
@@ -111,7 +113,7 @@ See also: `Model_CC`, `model_CC_Constructor()`, `reversibility()`
 
 """
 
-function swiftCC(ModelObject_CC::Model_CC, Tolerance::Float64=1e-6, printLevel::Int=1)
+function swiftCC(ModelObject_CC::Model_CC, Tolerance::Float64=1e-6, OctuplePrecision::Bool=false, printLevel::Int=1)
 
     ## Extract relevant information from the input model object
 
@@ -127,6 +129,7 @@ function swiftCC(ModelObject_CC::Model_CC, Tolerance::Float64=1e-6, printLevel::
 
     # Create an array of reaction IDs:
     Reaction_Ids = collect(1:n)
+
     irreversible_reactions_id, reversible_reactions_id = reversibility(lb, Reaction_Ids, 0)
     n_irr = length(irreversible_reactions_id)
     n_rev = length(reversible_reactions_id)
@@ -138,7 +141,13 @@ function swiftCC(ModelObject_CC::Model_CC, Tolerance::Float64=1e-6, printLevel::
     ub_u = ones(n_irr)
 
     # Create a new optimization model using the GLPK optimizer:
-    model = Model(GLPK.Optimizer)
+    if OctuplePrecision
+        model = GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat})
+        settings = Clarabel.Settings()
+        settings = Clarabel.Settings(verbose = false, time_limit = 5)
+    else
+        model = Model(GLPK.Optimizer)
+    end
 
     # Define variables V and u with their lower and upper bounds:
     @variable(model, lb[i] <= V[i = 1:n] <= ub[i])

src/ConsistencyChecking/TheNaiveApproach.jl

@@ -10,7 +10,9 @@ module TheNaiveApproach
 
 export find_blocked_reactions
 
-using GLPK, JuMP, COBREXA, Distributed
+using GLPK, JuMP, COBREXA, Distributed, Clarabel
+
+import CDDLib
 
 include("../Pre_Processing/Pre_processing.jl")
 
@@ -52,7 +54,7 @@ See also: `dataOfModel()`, `reversibility()`
 
 """
 
-function find_blocked_reactions(model::CoreModel, Tolerance::Float64=1e-6, printLevel::Int=1)
+function find_blocked_reactions(model::CoreModel, Tolerance::Float64=1e-6, OctuplePrecision::Bool=false, printLevel::Int=1)
 
     ## Export data from model
 
@@ -71,7 +73,7 @@ function find_blocked_reactions(model::CoreModel, Tolerance::Float64=1e-6, print
     printstyled("Homogenization:\n"; color=:cyan)
 
     # Set the maximum value for M:
-    M = 1000000.0
+    M = getM()
 
     ## Loop through each value in the array "lb" and "ub"
 
@@ -101,7 +103,14 @@ function find_blocked_reactions(model::CoreModel, Tolerance::Float64=1e-6, print
     irreversible_unblocked_reactions_id = []
 
     # Create a new optimization model using the GLPK optimizer:
-    model_irr = Model(GLPK.Optimizer)
+
+    if OctuplePrecision
+        model_irr = GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat})
+        settings = Clarabel.Settings()
+        settings = Clarabel.Settings(verbose = false, time_limit = 5)
+    else
+        model_irr = Model(GLPK.Optimizer)
+    end
 
     # Define the variable V for each reaction, with its lower and upper bounds:
     @variable(model_irr, lb[i] <= V[i = 1:n] <= ub[i])
@@ -143,7 +152,14 @@ function find_blocked_reactions(model::CoreModel, Tolerance::Float64=1e-6, print
     reversible_blocked_reactions_id = []
 
     # Create a new optimization model using the GLPK optimizer:
-    model_rev = Model(GLPK.Optimizer)
+
+    if OctuplePrecision
+        model_rev = GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat})
+        settings = Clarabel.Settings()
+        settings = Clarabel.Settings(verbose = false, time_limit = 5)
+    else
+        model_rev = Model(GLPK.Optimizer)
+    end
 
     # Define the variable V for each reaction, with its lower and upper bounds:
     @variable(model_rev, lb[i] <= V[i = 1:n] <= ub[i])

src/Pre_Processing/Pre_processing.jl

@@ -10,7 +10,9 @@ module Pre_processing
 
 export dataOfModel, getM, getTolerance, reversibility, check_duplicate_reactions, homogenization, remove_zeroRows, Model_Correction, model_Correction_Constructor, distributedReversibility_Correction
 
-using GLPK, JuMP, COBREXA, SparseArrays, Distributed, SharedArrays, Distributed
+using GLPK, JuMP, COBREXA, SparseArrays, Distributed, SharedArrays, Distributed, Clarabel
+
+import CDDLib
 
 """
     dataOfModel(model)
@@ -549,7 +551,7 @@ See also: `dataOfModel()`, `reversibility()`, `homogenization()`, 'getTolerance(
 
 """
 
-function distributedReversibility_Correction(ModelObject_Correction::Model_Correction, blocked_index_rev::Vector{Int64}, printLevel::Int=1)
+function distributedReversibility_Correction(ModelObject_Correction::Model_Correction, blocked_index_rev::Vector{Int64}, OctuplePrecision::Bool=false, printLevel::Int=1)
 
     ## Extract relevant information from the input model object
 
@@ -563,7 +565,7 @@ function distributedReversibility_Correction(ModelObject_Correction::Model_Corre
     n = length(lb)
 
     # Set the tolerance value:
-    Tolerance = 1e-6
+    Tolerance = getTolerance()
 
     # Initialize empty arrays to store the IDs of blocked reversible reactions in the forward and backward directions:
     rev_blocked_fwd = Array{Int64}([])
@@ -576,53 +578,65 @@ function distributedReversibility_Correction(ModelObject_Correction::Model_Corre
     # Iterate over all reversible reactions in the model:
     @sync @distributed for i in reversible_reactions_id
 
-        if i in blocked_index_rev
-            Correction[i,i] = +2.0
-            continue
+        # Create a local model object for each worker process:
+        if OctuplePrecision
+            local_model = GenericModel{BigFloat}(Clarabel.Optimizer{BigFloat})
+            settings = Clarabel.Settings()
+            settings = Clarabel.Settings(verbose = false, time_limit = 5)
         else
-
-            # Create a local model object for each worker process:
             local_model = Model(GLPK.Optimizer)
+        end
 
-            # Define variables V:
-            @variable(local_model, lb[i] <= V[i = 1:n] <= ub[i])
+        # Define variables V:
+        @variable(local_model, lb[i] <= V[i = 1:n] <= ub[i])
 
-            # Add the stoichiometric constraints to the model:
-            @constraint(local_model, S * V .== 0)
+        # Add the stoichiometric constraints to the model:
+        @constraint(local_model, S * V .== 0)
 
-            # Set the objective function to maximize the flux through reaction i in the forward direction:
-            @objective(local_model, Max, V[i])
+        # Set the objective function to maximize the flux through reaction i in the forward direction:
+        @objective(local_model, Max, V[i])
 
-            # Add a constraint that limits the flux through reaction i in the forward direction to be less than or equal to 1:
-            @constraint(local_model, V[i] <= 1)
+        # Add a constraint that limits the flux through reaction i in the forward direction to be less than or equal to 1:
+        @constraint(local_model, c_irr, V[i] <= 1)
 
-            # Optimize the model and retrieve the objective value:
-            optimize!(local_model)
-            opt_fwd = objective_value(local_model)
+        # Optimize the model and retrieve the objective value:
+        optimize!(local_model)
+        opt_fwd = objective_value(local_model)
 
-            # Set the objective function to minimize the flux through reaction i in the backward direction:
-            @objective(local_model, Min, V[i])
+        # Remove the current constraint from the model:
+        delete(local_model, c_irr)
 
-            # Add a constraint that limits the flux through reaction i in the backward direction to be greater than or equal to -1:
-            @constraint(local_model, V[i] >= -1)
+        # Unregister the current constraint from the model:
+        unregister(local_model, :c_irr)
 
-            # Optimize the model and retrieve the objective value:
-            optimize!(local_model)
-            opt_back = objective_value(local_model)
+        # Set the objective function to minimize the flux through reaction i in the backward direction:
+        @objective(local_model, Min, V[i])
 
-            # The reaction is considered to be blocked:
-            if isapprox(opt_fwd, 0, atol=Tolerance) && isapprox(opt_back, 0, atol=Tolerance)
-                Correction[i,i] = +2.0
-                continue
-            # If the objective value is approximately 0, the reaction is considered to be blocked in the forward direction:
-            elseif isapprox(opt_fwd, 0, atol=Tolerance)
-                Correction[i,i] = +1.0
-            # If the objective value is approximately 0, the reaction is considered to be blocked in the backward direction:
-            elseif isapprox(opt_back, 0, atol=Tolerance)
-                Correction[i,i] = -1.0
-            else
-                continue
-            end
+        # Add a constraint that limits the flux through reaction i in the backward direction to be greater than or equal to -1:
+        @constraint(local_model, c_rev, V[i] >= -1)
+
+        # Optimize the model and retrieve the objective value:
+        optimize!(local_model)
+        opt_back = objective_value(local_model)
+
+        # Remove the current constraint from the model:
+        delete(local_model, c_rev)
+
+        # Unregister the current constraint from the model:
+        unregister(local_model, :c_rev)
+
+        # The reaction is considered to be blocked:
+        if isapprox(opt_fwd, 0, atol=Tolerance) && isapprox(opt_back, 0, atol=Tolerance)
+            Correction[i,i] = +2.0
+            continue
+        # If the objective value is approximately 0, the reaction is considered to be blocked in the forward direction:
+        elseif isapprox(opt_fwd, 0, atol=Tolerance)
+            Correction[i,i] = +1.0
+        # If the objective value is approximately 0, the reaction is considered to be blocked in the backward direction:
+        elseif isapprox(opt_back, 0, atol=Tolerance)
+            Correction[i,i] = -1.0
+        else
+            continue
         end
 
     end

src/QFCA/distributedQFCA.jl

@@ -133,7 +133,7 @@ See also: `dataOfModel()`, `reversibility()`, `homogenization()`, `Model_QFCA`,
 
 """
 
-function distributedQFCA(ModelObject_QFCA::Model_QFCA, blocked_index::Vector{Int64}, removing::Bool=false, Tolerance::Float64=1e-6, printLevel::Int=1)
+function distributedQFCA(ModelObject_QFCA::Model_QFCA, blocked_index::Vector{Int64}, removing::Bool=false, Tolerance::Float64=1e-6, OctuplePrecision::Bool=false, printLevel::Int=1)
 
 
     ## Extract relevant information from the input model object
@@ -217,7 +217,7 @@ function distributedQFCA(ModelObject_QFCA::Model_QFCA, blocked_index::Vector{Int
 
             # Calculate the set of blocked reactions and dual variables for the modified network:
             model_CC_Constructor(ModelObject_CC ,S_noBlocked, Metabolites, Reactions_noBlocked, Genes, row_noBlocked, col_noBlocked, lb_noBlocked, ub_noBlocked)
-            blocked, ν = swiftCC(ModelObject_CC, Tolerance, 0)
+            blocked, ν = swiftCC(ModelObject_CC, Tolerance, OctuplePrecision, 0)
 
             # Update indices of blocked reactions after removing ith reaction:
             for j = 1:length(blocked)
@@ -248,7 +248,7 @@ function distributedQFCA(ModelObject_QFCA::Model_QFCA, blocked_index::Vector{Int
 
             # Find the set of blocked reactions and dual variables for the modified network:
             model_CC_Constructor(ModelObject_CC ,S_noBlocked, Metabolites, Reactions_noBlocked, Genes, row_noBlocked, col_noBlocked, lb_noBlocked, ub_noBlocked)
-            blocked, ν = swiftCC(ModelObject_CC, Tolerance, 0)
+            blocked, ν = swiftCC(ModelObject_CC, Tolerance, OctuplePrecision, 0)
 
             # Update DC_Matrix based on the blocked reactions:
             DC_Matrix[i,blocked] .= 1.0