added tests for hit and run sampler, broken parallel

src/analysis/sampling/hit_and_run.jl

@@ -7,7 +7,7 @@
         keepevery = _constants.sampling_keep_iters,
         samplesize = _constants.sampling_size,
         warmup_indices = collect(1:n_reactions(model)),
-        workers = [myid()],
+        workerids = [myid()],
         nchains = 1,
     )
 
@@ -60,7 +60,7 @@ chains = hit_and_run(
     N = 1000_000,
     nchains = 3,
     modifications = [change_constraint("EX_glc__D_e",-8, -8)]
-    workers = workers()
+    workerids = workers()
     )
 ```
 """
@@ -68,11 +68,11 @@ function hit_and_run(
     model,
     optimizer;
     modifications = [],
-    N = 1000,
+    N = 150_000,
     keepevery = _constants.sampling_keep_iters,
     samplesize = _constants.sampling_size,
     warmup_indices = collect(1:n_reactions(model)),
-    workers = [myid()],
+    workerids = [myid()],
     nchains = 1,
 )
 
@@ -81,16 +81,16 @@ function hit_and_run(
         model,
         optimizer;
         modifications = modifications,
-        workers = workers, # parallel
+        workerids = workerids, # parallel
         warmup_points = warmup_indices,
     )
 
     # load warmup points to workers
-    save_at.(workers, :cobrexa_ws, Ref(:($ws)))
-    save_at.(workers, :cobrexa_lbs, Ref(:($lbs)))
-    save_at.(workers, :cobrexa_ubs, Ref(:($ubs)))
+    save_at.(workerids, :cobrexa_ws, Ref(:($ws)))
+    save_at.(workerids, :cobrexa_lbs, Ref(:($lbs)))
+    save_at.(workerids, :cobrexa_ubs, Ref(:($ubs)))
 
-    # do in parallel! 
+    # sample in parallel! 
     samples = dpmap(
         x -> :($COBREXA._serial_hit_and_run(
             cobrexa_ws,
@@ -100,14 +100,14 @@ function hit_and_run(
             $keepevery,
             $N,
         )),
-        CachingPool(workers),
+        CachingPool(workerids),
         1:nchains,
     )
 
     # remove warmup points from workers
-    map(fetch, remove_from.(workers, :cobrexa_ws))
-    map(fetch, remove_from.(workers, :cobrexa_lbs))
-    map(fetch, remove_from.(workers, :cobrexa_ubs))
+    map(fetch, remove_from.(workerids, :cobrexa_ws))
+    map(fetch, remove_from.(workerids, :cobrexa_lbs))
+    map(fetch, remove_from.(workerids, :cobrexa_ubs))
 
     # not sure how to do this better - cat/vcat doesn't work, oh well 
     vals = zeros(samplesize, length(lbs), nchains)
@@ -159,7 +159,7 @@ function _serial_hit_and_run(ws, lbs, ubs, samplesize, keepevery, N)
         end
 
         if λmax <= λmin || λmin == -Inf || λmax == Inf # this sometimes can happen
-            #     @warn "Infeasible direction at iteration $(n)..." # noisy
+            #     @warn "Infeasible direction at iteration $(n)..." # TODO: make this optional/ add to a logger... very noisy otherwise
             continue
         end
 

test/analysis/sampling/hit_and_run.jl

@@ -7,22 +7,46 @@
 
     model = load_model(StandardModel, model_path)
 
+    using Distributed
+
+    nps = 4 - nprocs()
+
+    # load extra processes, have at least 4 available
+    if 1 <= nps <= 3 
+        addprocs(nps)
+    end
+    @everywhere using COBREXA, Tulip
+
+    # Serial test
     chains = hit_and_run(
         model,
         Tulip.Optimizer;
-        N = 5000_000,
+        N = 10_000,
         nchains = 3,
-        samplesize = 10_000,
+        samplesize = 1000,
         modifications = [
             change_constraint("EX_glc__D_e", -10, -10),
             change_constraint("BIOMASS_Ecoli_core_w_GAM", 0.8, 0.8),
         ],
+        workerids = [myid()],
     )
 
-    # # The sampling converges very slowly, so can't really do an accurate test 
-    # with so few samples
-    # this test is ugly, there must be a better way to get the mean
-    @test isapprox(mean(chains[[:PFL]]).nt.mean[1], 0.8, atol = 0.5)
+    # Do not test for convergence, that requires too many samples
+    @test isapprox(mean(chains[[:PFL]]).nt.mean[1], 0.8, atol = 1.0)
+
+    # parallel tests (still broken, not sure why)
+    chainsparallel = hit_and_run(
+        model,
+        Tulip.Optimizer;
+        N = 10_000,
+        nchains = 3,
+        samplesize = 1000,
+        modifications = [
+            change_constraint("EX_glc__D_e", -10, -10),
+            change_constraint("BIOMASS_Ecoli_core_w_GAM", 0.8, 0.8),
+        ],
+        workerids = workers(),
+    )
 
-    # TODO: parallel tests (still broken, not sure why)
+    @test isapprox(mean(chainsparallel[[:PFL]]).nt.mean[1], 0.8, atol = 1.0)
 end

test/analysis/sampling/warmup.jl

@@ -32,7 +32,7 @@
     @test ubs[ind] ≈ 4
     
     # Parallel test
-    ws, lbs, ubs = warmup(
+    wsparallel, lbsparallel, ubsparallel = warmup(
         model,
         Tulip.Optimizer;
         modifications = [change_constraint("EX_glc__D_e", -4, 4)],
@@ -41,8 +41,8 @@
     )
 
     ind = first(indexin(["EX_glc__D_e"], reactions(model)))
-    @test size(ws) == (95, 2)
-    @test size(ws[1,1]) == (95,)
-    @test lbs[ind] ≈ -4
-    @test ubs[ind] ≈ 4
+    @test size(wsparallel) == (95, 2)
+    @test size(wsparallel[1,1]) == (95,)
+    @test lbsparallel[ind] ≈ -4
+    @test ubsparallel[ind] ≈ 4
 end