Merge pull request #33 from iManGHD/master

Coupling Coefficients

Project.toml

@@ -13,6 +13,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
 
 [compat]
 COBREXA = "1.0, 1.2"

example/sparseQFCA.ipynb

@@ -2,23 +2,15 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "┌ Info: Recompiling stale cache file C:\\Users\\mtefa\\.julia\\compiled\\v1.0\\sparseQFCA\\Q0hmV.ji for sparseQFCA [47e41d82-2d4b-11e9-01da-5b8950e38350]\n",
-      "└ @ Base loading.jl:1190\n"
-     ]
-    },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "  2.698186 seconds (29.77 M allocations: 790.907 MiB, 11.91% gc time)\n",
-      "961.782199 seconds (3.19 G allocations: 411.142 GiB, 12.86% gc time)\n",
+      "  2.841473 seconds (29.15 M allocations: 760.187 MiB, 13.26% gc time, 23.18% compilation time)\n",
+      "861.983010 seconds (3.91 G allocations: 351.983 GiB, 7.58% gc time, 1.66% compilation time)\n",
       "The answer is correct.\n",
       "The answer is correct.\n"
      ]
@@ -59,15 +51,15 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.0.3",
+   "display_name": "Julia 1.7.2",
    "language": "julia",
-   "name": "julia-1.0"
+   "name": "julia-1.7"
   },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.0.3"
+   "version": "1.7.2"
   }
  },
  "nbformat": 4,

src/Consistency Checking/SwiftCC.jl

@@ -1,55 +1,110 @@
 #-------------------------------------------------------------------------------------------
-
 #=
     Purpose:    Finding blocked reactions in metabolic network
     Author:     Iman Ghadimi, Mojtaba Tefagh - Sharif University of Technology - Iran
     Date:       July 2022
 =#
-
 #-------------------------------------------------------------------------------------------
 
 module SwiftCC
-export swiftCC
-
-using Distributed
+export MyModel, myModel_Constructor, swiftCC
 
-@everywhere using GLPK, JuMP, COBREXA, LinearAlgebra, SparseArrays
+using GLPK, JuMP, COBREXA, LinearAlgebra, SparseArrays, Distributed
 
 include("../Data Processing/pre_processing.jl")
 
 using .pre_processing
 
+
+"""
+    MyModel(S, Metabolites, Reactions, Genes, m, n, lb, ub)
+
+A general type for storing a StandardModel which contains the following fields:
+
+- `S`:              LHS matrix (m x n)
+- `Metabolites`:    List of metabolic network metabolites.
+- `Reactions`:      List of metabolic network reactions.
+- `Genes`:          List of metabolic network reactions.
+- `m`:              Number of rows of stoichiometric matrix.
+- `n`:              Number of columns of stoichiometric matrix.
+- `lb`:             Lower bound vector (n x 1)
+- `ub`:             Upper bound vector (n x 1)
+
+"""
+
+mutable struct MyModel
+    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}
+end
+
+"""
+    myModel_Constructor(ModelObject, S, Metabolites, Reactions, Genes, m, n, lb, ub)
+
+A function that initializes a newly created object of MyModel.
+
+# INPUTS
+
+-'ModelObject'      a newly object of MyModel.
+- `S`:              LHS matrix (m x n)
+- `Metabolites`:    List of metabolic network metabolites.
+- `Reactions`:      List of metabolic network reactions.
+- `Genes`:          List of metabolic network reactions.
+- `m`:              Number of rows of stoichiometric matrix.
+- `n`:              Number of columns of stoichiometric matrix.
+- `lb`:             Lower bound vector (n x 1)
+- `ub`:             Upper bound vector (n x 1)
+
+"""
+
+function myModel_Constructor(ModelObject::MyModel, 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})
+     ModelObject.S = S
+     ModelObject.Metabolites = Metabolites
+     ModelObject.Reactions = Reactions
+     ModelObject.Genes = Genes
+     ModelObject.m = m
+     ModelObject.n = n
+     ModelObject.lb = lb
+     ModelObject.ub = ub
+end
 """
     swiftCC(ModelObject)
 
 A function that finds blocked reactions for a metabolic network.
 
 # INPUTS
 
-- `ModelObject`:        a newly object of MyModel.
+- `ModelObject`:        A newly object of MyModel.
 
 # OPTIONAL INPUTS
 
--
+- `Tolerance`:          A small number that represents the level of error tolerance.
 
 # OUTPUTS
 
-- `blocked_names`:      Blocked reaction Names.
+- `blocked_index`:      Index of blocked reactions.
+- `dualVar`:            Dual variables of a specific constraint.
 
 # EXAMPLES
 
 - Full input/output example
 ```julia
-julia> blocked_reactions = swiftCC(ModelObject)
+julia> blocked_index, dual_var = swiftCC(ModelObject)
 ```
 
-See also: `MyModel`, myModel_Constructor(), 'getTolerance()', `reversibility()`, `homogenization()`
+See also: `MyModel`, myModel_Constructor(), `reversibility()`, `homogenization()`
 
 """
 
-@everywhere function swiftCC(ModelObject::MyModel)
+function swiftCC(ModelObject::MyModel, Tolerance::Float64=1e-6)
 
-    # Exporting data from ModelObject:
+    ## Export data from ModelObject:
 
     S = ModelObject.S
     Metabolites = ModelObject.Metabolites
@@ -59,97 +114,104 @@ See also: `MyModel`, myModel_Constructor(), 'getTolerance()', `reversibility()`,
     lb = ModelObject.lb
     ub = ModelObject.ub
 
-    # assigning a small value to Tolerance representing the level of error tolerance:
-
-    Tolerance = getTolerance()
-
-    # Determining the reversibility of a reaction:
+    ## Determine the reversibility of a reaction:
 
     irreversible_reactions_id, reversible_reactions_id = reversibility(lb)
 
-    # Determining the number of irreversible and reversible reactions:
+    ## Determine the number of irreversible and reversible reactions:
 
     n_irr = length(irreversible_reactions_id)
     n_rev = length(reversible_reactions_id)
 
-    # Calculating the dimensions of the S matrix:
+    ## Calculate the dimensions of the S matrix:
 
     row_num, col_num = size(S)
 
-    # Finding irreversible blocked reactions in 1 LP problem:
+    ## Find irreversible blocked reactions in 1 LP problem:
 
+    # Set the lower and upper bounds of the "u" variables to 0 and 1, respectively
     lb_u = zeros(n_irr)
     ub_u = ones(n_irr)
+
+    # Create a new optimization model using the GLPK optimizer
     model = Model(GLPK.Optimizer)
+
+    # Define variables V and u with their lower and upper bounds
     @variable(model, lb[i] <= V[i = 1:n] <= ub[i])
     @variable(model, lb_u[i] <= u[i = 1:n_irr] <= ub_u[i])
-    @constraint(model, S * V .== 0)
+
+    # Add a constraint that ensures the mass balance of the system
+    @constraint(model, c1, S * V .== 0)
+
+    # Define the objective function to be the sum of all "u" variables
     objective_function = ones(n_irr)'*u
+
+    # Set the optimization objective to maximize the objective function
     @objective(model, Max, objective_function)
+
+    # Add a constraint for each irreversible reaction that sets its "u" variable to be less than or equal to its corresponding V variable
     @constraint(model, [i in 1:n_irr], u[i] <= V[irreversible_reactions_id[i]])
+
+    # Solve the optimization problem
     optimize!(model)
 
+    # Get the dual variables for the stoichiometric constraint
+    dualVar = dual.(c1)
+
+    # Initialize an empty list to store the IDs of blocked irreversible reactions
     irr_blocked_reactions = []
+
+    # Iterate over the irreversible reactions and check if the corresponding u variable is close to 0, If it is, the reaction is considered blocked and its ID is added to the list
     for i in range(1,n_irr; step=1)
         if isapprox(value(u[i]), 0.0, atol = Tolerance)
             append!(irr_blocked_reactions, irreversible_reactions_id[i])
         end
     end
 
-    # Determining the number of irreversible blocked reactions:
+    # Determine the number of irreversible blocked reactions
 
     irr_blocked_num = length(irr_blocked_reactions)
 
-    # Finding reversible blocked reactions:
-
-    # S_Transpose:
+    ## Find reversible blocked reactions:
 
+    # Creating S_transpose
     S_transpose = S'
 
-    # Calculating the dimensions of the S_transpose matrix:
-
+    # Determining the dimensions of the S_transpose matrix
     row_num_trans, col_num_trans = size(S_transpose)
 
-    # Removing irrevesible blocked from Stoichiometric Matrix:
-
+    # Removing irreversibly blocked reactions from the Stoichiometric Matrix
     S_transpose_noIrrBlocked = S_transpose[setdiff(1:end, irr_blocked_reactions), :]
 
-    # Constructing the I_reversible Matrix:
-
+    # Creating the I_reversible Matrix:
     I_reversible = zeros(n, n_rev)
 
+    # Populating the I_reversible Matrix with 1 in the rows corresponding to reversible reactions
     rev_id = 1
     for col in eachcol(I_reversible)
         col[reversible_reactions_id[rev_id]] = 1.0
-        rev_id = rev_id + 1
+        rev_id += 1
     end
 
-    # Removing irrevesible blocked from I_reversible Matrix:
-
+    # Removing irreversibly blocked reactions from the I_reversible Matrix
     I_reversible = I_reversible[setdiff(1:end, irr_blocked_reactions), :]
 
-    # Calculating the dimensions of the S_transpose_noIrrBlocked and I_reversible matrix :
-
+    # Determining the dimensions of the S_transpose_noIrrBlocked and I_reversible matrices
     S_trans_row, S_trans_col = size(S_transpose_noIrrBlocked)
     I_row, I_col = size(I_reversible)
 
-    # using Gaussian elimination to find reversible blocked reactions:
-
+    # Solving the system of equations using Gaussian elimination to identify blocked reversible reactions
     X = S_transpose_noIrrBlocked \ I_reversible
-
     Sol = (S_transpose_noIrrBlocked*X) - I_reversible
 
-    # Calculating the dimensions of the S_transpose_noIrrBlocked matrix:
-
+    # Determining the dimensions of the Sol matrix
     row_sol, col_sol = size(Sol)
 
-    # Finding specific columns in Sol:
-
+    # Finding columns in Sol that correspond to blocked reversible reactions
     c = 0
-
     rev_blocked_reactions_col = []
     for col in eachcol(Sol)
-        c = c + 1
+        c += 1
         if isapprox(norm(col), 0, atol = Tolerance)
             append!(rev_blocked_reactions_col, c)
         end
@@ -160,15 +222,15 @@ See also: `MyModel`, myModel_Constructor(), 'getTolerance()', `reversibility()`,
         append!(rev_blocked_reactions, reversible_reactions_id[i])
     end
 
-    # Merging reversbile blocked reactions and irreversible blocked reactions:
+    # Union reversbile blocked reactions and irreversible blocked reactions
 
     blocked_index = []
     blocked_index = union(rev_blocked_reactions, irr_blocked_reactions)
 
-    # Returning a list consist of the Ids of the blocked reactions:
+    # Returning a list consist of the Ids of the blocked reactions
 
     blocked_index = sort(blocked_index)
-    return blocked_index
+    return blocked_index, dualVar
 end
 
 end