From 40ab1f9212bbd777f3a8b39ad8ef998672a143a2 Mon Sep 17 00:00:00 2001
From: stormaths <hackeint@gmail.com>
Date: Mon, 24 Jul 2023 16:54:54 +0100
Subject: [PATCH] Modified the parameters used to call the Gaussian and
 exponential HMC walks in bindings. Added tests to check they work as expected

---
 dingo/bindings/bindings.cpp     |  7 ++-
 tests/sampling_no_multiphase.py | 75 +++++++++++++++++++++++++++++++++
 2 files changed, 78 insertions(+), 4 deletions(-)
 create mode 100644 tests/sampling_no_multiphase.py

diff --git a/dingo/bindings/bindings.cpp b/dingo/bindings/bindings.cpp
index d4d70fbe..d3ec9292 100644
--- a/dingo/bindings/bindings.cpp
+++ b/dingo/bindings/bindings.cpp
@@ -104,10 +104,6 @@ double HPolytopeCPP::apply_sampling(int walk_len,
    starting_point = inner_point2;
    std::list<Point> rand_points;
    NT variance = 1.0;
-   NT a = NT(1)/(NT(2)*variance);
-   int dim = 50;
-   Point c(dim);
-   c = GetDirection<Point>::apply(dim, rng, false);
 
    if (method == 1) { // cdhr
       uniform_sampling<CDHRWalk>(rand_points, HP, rng, walk_len, number_of_points,
@@ -133,10 +129,13 @@ double HPolytopeCPP::apply_sampling(int walk_len,
                                    starting_point, number_of_points_to_burn);
    } 
    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);
    } 
diff --git a/tests/sampling_no_multiphase.py b/tests/sampling_no_multiphase.py
new file mode 100644
index 00000000..f46f6906
--- /dev/null
+++ b/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()