add parFVA

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 Clp = "e2554f3b-3117-50c0-817c-e040a3ddf72d"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"

src/COBREXA.jl

@@ -6,6 +6,7 @@ using DelimitedFiles
 using LinearAlgebra
 using JuMP
 using MAT
+using Distributed
 import Pkg
 
 # import src files
@@ -43,7 +44,7 @@ export speye, LinearModel, nReactions, nMetabolites, nCouplingConstraints,
        removeCouplingConstraints!, removeCouplingConstraints,
        changeCouplingBounds!, changeCouplingBounds,
        verifyConsistency, findExchangeReactions, findExchangeMetabolites,
-       solveLP, loadModel, fluxBalanceAnalysis, fluxVariabilityAnalysis,
-       writeModel, convertToExportable
+       solveLP, loadModel, fluxBalanceAnalysis, fluxVariabilityAnalysis, parFVA,
+       writeModel, convertToExportable, createParPool
 
 end

src/analysis/fva.jl

@@ -39,3 +39,51 @@ function fluxVariabilityAnalysis(model::LinearModel, reactions::Array{Int64, 1},
    end
    return fluxes
 end
+
+"""
+(Local) multi-process version of FVA that assigns the reactions to several
+processes running in parallel
+`workers` should be a list of process ids returned by `createParPool`
+"""
+function parFVA(model::LinearModel, reactions::Array{Int64, 1}, solverName::Symbol, workersToUse::Array{Int64,1})
+   nReacs = length(reactions)
+   nWorkers = length(workersToUse)
+   if nReacs < nWorkers
+      @info "Number of workers exceeds number of reactions. 1 worker per reaction will be used."
+      workersToUse = workersToUse[1:nReacs]
+      nWorkers = nReacs
+   end
+   remRefs = Array{Future}(undef, nWorkers)
+   fluxes = zeros(nReacs, 2)
+   alloReacs = allocateReacs(reactions, nWorkers)
+
+   @sync for (round, pid) in enumerate(workersToUse)
+      @async remRefs[round] = @spawnat pid begin
+         fluxVariabilityAnalysis(model, alloReacs[round], eval(solverName).Optimizer)
+      end
+   end
+   @sync for (round, pid) in enumerate(workersToUse)
+      fluxes[alloReacs[round], :] = fetch(remRefs[round])
+   end
+
+   return fluxes
+end
+
+"""
+Auxiliary function for `parFVA` to divide the list of reactions evenly (for now)
+"""
+function allocateReacs(reactions::Array{Int64, 1}, nWorkers::Int)
+   nReacs = length(reactions)
+   steps = floor(Int, nReacs/nWorkers) * ones(Int, nWorkers)
+   steps[1:nReacs%nWorkers] = steps[1:nReacs%nWorkers] .+ 1
+   steps = cumsum(steps)
+   allocatedReacs = Vector(undef, nWorkers)
+
+   startI = 1
+   for (i, endI) in enumerate(steps)
+      allocatedReacs[i] = reactions[startI:endI]
+      startI = endI + 1
+   end
+
+   return allocatedReacs
+end

src/base/utilities.jl

@@ -15,3 +15,12 @@ function Base.copy(model::LinearModel)
     return LinearModel(model.S, model.b, model.C, model.cl, model.cu,
                        model.c, model.xl, model.xu, model.rxns, model.mets)
 end
+
+"""
+Create a pool of parallel workers to pass to, eg, `parFVA`
+"""
+function createParPool(nWorkers)
+    addprocs(nWorkers, topology=:master_worker)
+
+    return workers()[end-nWorkers+1: end]
+end

test/analysis/fva.jl

@@ -12,3 +12,12 @@
     @test size(fluxes) == (1, 2)
     @test fluxes == Array{Float64, 2}([2 2])
 end
+
+@testset "Parallel FVA" begin
+    cp = test_simpleLP()
+    pids = createParPool(2)
+    @everywhere using COBREXA, GLPK
+    fluxes = parFVA(cp, [1;2], :GLPK, pids)
+    @test fluxes ≈ [1. 1.;
+                    2. 2.]
+end

test/base/utilities.jl

@@ -13,3 +13,10 @@
     cp = test_coupledLP()
     @test nCouplingConstraints(cp) == 2000
 end
+
+@testset "Allocating reactions" begin
+    @test COBREXA.allocateReacs([1], 1) == [[1]]
+    @test COBREXA.allocateReacs(collect(1:3), 2) == [[1;2], [3]]
+    @test COBREXA.allocateReacs(collect(1:8), 4) == [[1;2], [3;4], [5;6], [7;8]]
+    @test COBREXA.allocateReacs(collect(1:11), 3) == [[1;2;3;4], [5;6;7;8], [9;10;11]]
+end

test/runtests.jl

@@ -6,6 +6,7 @@ using GLPK
 using Clp
 using COBREXA
 using MAT
+using Distributed
 
 """
     runSuite(baseDir)