Add mmcs method to sampling (#108)

* Upgrade volesti and pass method as string avoiding uneeded integer checks

* Add mmcs method to sampling

* Add tests that call mmcs from sampling function

* Update volesti submodule

---------

Co-authored-by: akisschinas <[email protected]>

Author: akisschinas

build.py

@@ -23,7 +23,7 @@ def build(setup_kwargs):
 
         # gcc arguments hack: enable optimizations
         os.environ["CFLAGS"] = [
-            "-std=c++11",
+            "-std=c++17",
             "-O3",
             "-DBOOST_NO_AUTO_PTR",
             "-ldl",

dingo/PolytopeSampler.py

@@ -159,9 +159,9 @@ def generate_steady_states(
         self._T_shift = np.add(self._T_shift, Tr_shift)
 
         return steady_states
-    
+
     def generate_steady_states_no_multiphase(
-        self, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None
+        self, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, ess=1000
     ):
         """A member function to sample steady states.
 
@@ -171,17 +171,17 @@ def generate_steady_states_no_multiphase(
         burn_in -- the number of points to burn before sampling
         thinning -- the walk length of the chain
         """
-        	
+
         self.get_polytope()
 
         P = HPolytope(self._A, self._b)
-        
+
         if bias_vector is None:
             bias_vector = np.ones(self._A.shape[1], dtype=np.float64)
         else:
             bias_vector = bias_vector.astype('float64')
 
-        samples = P.generate_samples(method, n, burn_in, thinning, variance, bias_vector, self._parameters["solver"])
+        samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, self._parameters["solver"], ess)
         samples_T = samples.T
 
         steady_states = map_samples_to_steady_states(
@@ -213,10 +213,10 @@ def sample_from_polytope(
 
 
         return samples
-    
+
     @staticmethod
     def sample_from_polytope_no_multiphase(
-        A, b, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, solver=None
+        A, b, method = 'billiard_walk', n=1000, burn_in=0, thinning=1, variance=1.0, bias_vector=None, solver=None, ess=1000
     ):
         """A static function to sample from a full dimensional polytope with an MCMC method.
 
@@ -232,10 +232,10 @@ def sample_from_polytope_no_multiphase(
             bias_vector = np.ones(A.shape[1], dtype=np.float64)
         else:
             bias_vector = bias_vector.astype('float64')
-            
+
         P = HPolytope(A, b)
 
-        samples = P.generate_samples(method, n, burn_in, thinning, variance, bias_vector, solver)
+        samples = P.generate_samples(method.encode('utf-8'), n, burn_in, thinning, variance, bias_vector, solver, ess)
 
         samples_T = samples.T
         return samples_T
@@ -246,7 +246,7 @@ def round_polytope(
     ):
         P = HPolytope(A, b)
         A, b, Tr, Tr_shift, round_value = P.rounding(method, solver)
-        
+
         return A, b, Tr, Tr_shift
 
     @staticmethod

dingo/bindings/bindings.cpp

@@ -38,7 +38,7 @@
 // for rounding
 #include "preprocess/min_sampling_covering_ellipsoid_rounding.hpp"
 #include "preprocess/svd_rounding.hpp"
-#include "preprocess/max_inscribed_ellipsoid_rounding.hpp"
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
 
 typedef double NT;
 typedef Cartesian<NT>    Kernel;
@@ -140,7 +140,8 @@ class HPolytopeCPP{
 
       // the apply_sampling() function
       double apply_sampling(int walk_len, int number_of_points, int number_of_points_to_burn,
-                            int method, double* inner_point, double radius, double* samples, double variance_value, double* bias_vector);
+                            char* method, double* inner_point, double radius, double* samples,
+                            double variance_value, double* bias_vector, int ess);
 
       void mmcs_initialize(int d, int ess, bool psrf_check, bool parallelism, int num_threads);
 

dingo/bindings/bindings.h

@@ -1,6 +1,6 @@
 # This is a cython wrapper for the C++ library volesti
 # volesti (volume computation and sampling library)
-  
+
 # Copyright (c) 2012-2021 Vissarion Fisikopoulos
 # Copyright (c) 2018-2021 Apostolos Chalkis
 # Copyright (c) 2020-2021 Pedro Zuidberg Dos Martires
@@ -54,8 +54,9 @@ cdef extern from "bindings.h":
 
       # Random sampling
       double apply_sampling(int walk_len, int number_of_points, int number_of_points_to_burn, \
-                              int method, double* inner_point, double radius, double* samples, double variance_value, double* bias_vector)
-      
+                            char* method, double* inner_point, double radius, double* samples, \
+                            double variance_value, double* bias_vector, int ess)
+
       # Initialize the parameters for the (m)ultiphase (m)onte (c)arlo (s)ampling algorithm
       void mmcs_initialize(unsigned int d, int ess, int psrf_check, int parallelism, int num_threads);
 
@@ -117,45 +118,22 @@ cdef class HPolytope:
          raise Exception('"{}" is not implemented to compute volume. Available methods are: {}'.format(vol_method, volume_methods))
 
    # Likewise, the generate_samples() function
-   def generate_samples(self, method, number_of_points, number_of_points_to_burn, walk_len, variance_value, bias_vector, solver = None):
+   def generate_samples(self, method, number_of_points, number_of_points_to_burn, walk_len,
+                        variance_value, bias_vector, solver = None, ess = 1000):
 
       n_variables = self._A.shape[1]
       cdef double[:,::1] samples = np.zeros((number_of_points, n_variables), dtype = np.float64, order = "C")
 
       # Get max inscribed ball for the initial polytope
       temp_center, radius = inner_ball(self._A, self._b, solver)
-      
+
       cdef double[::1] inner_point_for_c = np.asarray(temp_center)
-      
+
       cdef double[::1] bias_vector_ = np.asarray(bias_vector)
-      
-      if method == 'cdhr':
-         int_method = 1
-      elif method == 'rdhr':
-         int_method = 2
-      elif method == 'billiard_walk':
-         int_method = 3
-      elif method == 'ball_walk':
-         int_method = 4
-      elif method == 'dikin_walk':
-         int_method = 5
-      elif method == 'john_walk':
-         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
-      elif method == 'hmc_leapfrog_gaussian':
-         int_method = 10
-      elif method == 'hmc_leapfrog_exponential':
-         int_method = 11
-      else:
-         raise RuntimeError("Uknown MCMC sampling method")
-      
+
       self.polytope_cpp.apply_sampling(walk_len, number_of_points, number_of_points_to_burn, \
-                                       int_method, &inner_point_for_c[0], radius, &samples[0,0], variance_value, &bias_vector_[0])
+                                       method, &inner_point_for_c[0], radius, &samples[0,0], \
+                                       variance_value, &bias_vector_[0], ess)
       return np.asarray(samples)
 
 
@@ -172,12 +150,12 @@ cdef class HPolytope:
       cdef double[:,::1] T_matrix = np.zeros((n_variables, n_variables), dtype=np.float64, order="C")
       cdef double[::1] shift = np.zeros((n_variables), dtype=np.float64, order="C")
       cdef double round_value
-      
+
       # Get max inscribed ball for the initial polytope
       center, radius = inner_ball(self._A, self._b, solver)
-      
+
       cdef double[::1] inner_point_for_c = np.asarray(center)
-      
+
       if rounding_method == 'john_position':
          int_method = 1
       elif rounding_method == 'isotropic_position':
@@ -190,8 +168,8 @@ cdef class HPolytope:
       self.polytope_cpp.apply_rounding(int_method, &new_A[0,0], &new_b[0], &T_matrix[0,0], &shift[0], round_value, &inner_point_for_c[0], radius)
 
       return np.asarray(new_A),np.asarray(new_b),np.asarray(T_matrix),np.asarray(shift),np.asarray(round_value)
-   
-   
+
+
    # (m)ultiphase (m)onte (c)arlo (s)ampling algorithm to generate steady states of a metabolic network
    def mmcs(self, ess = 1000, psrf_check = True, parallelism = False, num_threads = 2, solver = None):
 
@@ -205,7 +183,7 @@ cdef class HPolytope:
       cdef int N_ess = ess
       cdef bint check_psrf = bool(psrf_check) # restrict variables to {0,1} using Python's rules
       cdef bint parallel = bool(parallelism)
-      
+
       self.polytope_cpp.mmcs_initialize(n_variables, ess, check_psrf, parallel, num_threads)
 
       # Get max inscribed ball for the initial polytope
@@ -215,14 +193,14 @@ cdef class HPolytope:
       while True:
 
          check = self.polytope_cpp.mmcs_step(&inner_point_for_c[0], radius, N_samples)
-         
+
          if check > 1.0 and check < 2.0:
             break
 
          self.polytope_cpp.get_polytope_as_matrices(&new_A[0,0], &new_b[0])
          new_temp_c, radius = inner_ball(np.asarray(new_A), np.asarray(new_b), solver)
          inner_point_for_c = np.asarray(new_temp_c)
-      
+
       cdef double[:,::1] samples = np.zeros((n_variables, N_samples), dtype=np.float64, order="C")
       self.polytope_cpp.get_mmcs_samples(&T_matrix[0,0], &T_shift[0], &samples[0,0])
       self.polytope_cpp.get_polytope_as_matrices(&new_A[0,0], &new_b[0])

dingo/volestipy.pyx

@@ -27,7 +27,7 @@
 
 source_directory_list = ["dingo", join("dingo", "bindings")]
 
-compiler_args = ["-std=c++11", "-O3", "-DBOOST_NO_AUTO_PTR", "-ldl", "-lm", "-fopenmp"]
+compiler_args = ["-std=c++17", "-O3", "-DBOOST_NO_AUTO_PTR", "-ldl", "-lm", "-fopenmp"]
 lp_solve_compiler_args = ["-DYY_NEVER_INTERACTIVE", "-DLoadInverseLib=0", "-DLoadLanguageLib=0",
 "-DRoleIsExternalInvEngine", "-DINVERSE_ACTIVE=3", "-DLoadableBlasLib=0"]
 

setup.py

@@ -14,108 +14,61 @@
 from dingo import MetabolicNetwork, PolytopeSampler
 from dingo.pyoptinterface_based_impl import set_default_solver
 
+def sampling(model, testing_class):
+    sampler = PolytopeSampler(model)
+
+    #Gaussian hmc sampling
+    steady_states = sampler.generate_steady_states_no_multiphase(method = 'mmcs', ess=1000)
+
+    testing_class.assertTrue( steady_states.shape[0] == 95 )
+    testing_class.assertTrue( steady_states.shape[1] == 1000 )
+
+    #Gaussian hmc sampling
+    steady_states = sampler.generate_steady_states_no_multiphase(method = 'gaussian_hmc_walk', n=500)
+
+    testing_class.assertTrue( steady_states.shape[0] == 95 )
+    testing_class.assertTrue( steady_states.shape[1] == 500 )
+
+    #exponential hmc sampling
+    steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=500, variance=50)
+
+    testing_class.assertTrue( steady_states.shape[0] == 95 )
+    testing_class.assertTrue( steady_states.shape[1] == 500 )
+
+    #hmc sampling with Gaussian distribution
+    steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=500)
+
+    testing_class.assertTrue( steady_states.shape[0] == 95 )
+    testing_class.assertTrue( steady_states.shape[1] == 500 )
+
+    #hmc sampling with exponential distribution
+    steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=500, variance=50)
+
+    testing_class.assertTrue( steady_states.shape[0] == 95 )
+    testing_class.assertTrue( steady_states.shape[1] == 500 )
+
+    #steady_states[12].mean() seems to have a lot of discrepancy between experiments, so we won't check the mean for now
+    #self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
 
 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', n=500)
-       
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #exponential hmc sampling
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with Gaussian distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=500)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with exponential distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #steady_states[12].mean() seems to have a lot of discrepancy between experiments, so we won't check the mean for now
-        #self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+        sampling(model, self)
 
     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', n=500)
-        
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #exponential hmc sampling
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with Gaussian distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=500)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with exponential distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #steady_states[12].mean() seems to have a lot of discrepancy between experiments, so we won't check the mean for now
-        #self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
-
+        sampling(model, self)
 
     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', n=500)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #exponential hmc sampling
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'exponential_hmc_walk', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with Gaussian distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_gaussian', n=500)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #hmc sampling with exponential distribution
-        steady_states = sampler.generate_steady_states_no_multiphase(method = 'hmc_leapfrog_exponential', n=500, variance=50)
-
-        self.assertTrue( steady_states.shape[0] == 95 )
-        self.assertTrue( steady_states.shape[1] == 500 )
-        
-        #steady_states[12].mean() seems to have a lot of discrepancy between experiments, so we won't check the mean for now
-        #self.assertTrue( abs( steady_states[12].mean()  - 2.504 ) < 1e-03 )
+        sampling(model, self)