Exposed two HMC walks with exponential and Gaussian sampling from Volesti

dingo/bindings/bindings.cpp

@@ -103,6 +103,7 @@ double HPolytopeCPP::apply_sampling(int walk_len,
    HP.set_InnerBall(CheBall);
    starting_point = inner_point2;
    std::list<Point> rand_points;
+   NT variance = 1.0;
 
    if (method == 1) { // cdhr
       uniform_sampling<CDHRWalk>(rand_points, HP, rng, walk_len, number_of_points,
@@ -120,13 +121,26 @@ double HPolytopeCPP::apply_sampling(int walk_len,
    } else if (method == 5) { // dikin walk {
       uniform_sampling<DikinWalk>(rand_points, HP, rng, walk_len, number_of_points,
                                   starting_point, number_of_points_to_burn);
-   } else if (method == 6) { // dikin walk {
+   } else if (method == 6) { // john walk {
       uniform_sampling<JohnWalk>(rand_points, HP, rng, walk_len, number_of_points,
                                  starting_point, number_of_points_to_burn);
-   } else if (method == 7) { // dikin walk {
+   } else if (method == 7) { // vaidya walk {
       uniform_sampling<VaidyaWalk>(rand_points, HP, rng, walk_len, number_of_points,
                                    starting_point, number_of_points_to_burn);
-   } else {
+   } 
+   else if (method == 8) { // gaussian sampling with gaussian HMC exact walk {
+   NT a = NT(1)/(NT(2)*variance);
+   gaussian_sampling<GaussianHamiltonianMonteCarloExactWalk>(rand_points, HP, rng, walk_len, number_of_points, a,
+                                   starting_point, number_of_points_to_burn);
+   } 
+   else if (method == 9) { // exponential sampling with exponential HMC exact walk {
+   Point c(d);
+   c = GetDirection<Point>::apply(d, rng, false);
+   exponential_sampling<ExponentialHamiltonianMonteCarloExactWalk>(rand_points, HP, rng, walk_len, number_of_points, c, variance,
+                                   starting_point, number_of_points_to_burn);
+   } 
+
+   else {
       throw std::runtime_error("This function must not be called.");
    }
 

dingo/volestipy.pyx

@@ -149,6 +149,10 @@ cdef class HPolytope:
          int_method = 6
       elif method == 'vaidya_walk':
          int_method = 7
+      elif method == 'gaussian_hmc_walk':
+         int_method = 8
+      elif method == 'exponential_hmc_walk':
+         int_method = 9
       else:
          raise RuntimeError("Uknown MCMC sampling method")
 

tests/sampling_no_multiphase.py

@@ -0,0 +1,75 @@
+# dingo : a python library for metabolic networks sampling and analysis
+# dingo is part of GeomScale project
+
+# Copyright (c) 2022 Apostolos Chalkis
+# Copyright (c) 2022 Vissarion Fisikopoulos
+# Copyright (c) 2022 Haris Zafeiropoulos
+
+# Licensed under GNU LGPL.3, see LICENCE file
+
+import unittest
+import os
+from dingo import MetabolicNetwork, PolytopeSampler
+
+
+class TestSampling(unittest.TestCase):
+
+    def test_sample_json(self):
+
+        input_file_json = os.getcwd() + "/ext_data/e_coli_core.json"
+        model = MetabolicNetwork.from_json( input_file_json )
+        sampler = PolytopeSampler(model)
+
+        #gaussian hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+        #exponential hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+    def test_sample_mat(self):
+
+        input_file_mat = os.getcwd() + "/ext_data/e_coli_core.mat"
+        model = MetabolicNetwork.from_mat(input_file_mat)
+        sampler = PolytopeSampler(model)
+
+        #gaussian hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+        #exponential hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+
+    def test_sample_sbml(self):
+
+        input_file_sbml = os.getcwd() + "/ext_data/e_coli_core.xml"
+        model = MetabolicNetwork.from_sbml( input_file_sbml )
+        sampler = PolytopeSampler(model)
+
+        #gaussian hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+        #exponential hmc sampling
+        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk')
+
+        self.assertTrue( steady_states.shape[0] == 95 )
+        self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+
+
+
+if __name__ == "__main__":
+    unittest.main()