HMC Reflections test

test/CMakeLists.txt

@@ -348,6 +348,12 @@ add_test(NAME test_new_rdhr_uniform_MT COMMAND matrix_sampling_test -tc=new_rdhr
 add_test(NAME test_new_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_billiard_uniform_MT)
 add_test(NAME test_new_accelerated_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_accelerated_billiard_uniform_MT)
 
+add_executable(gaussian_hmc_reflections_test gaussian_hmc_reflections_test.cpp $<TARGET_OBJECTS:test_main>)
+add_test(NAME ghmc_reflections_test_cube COMMAND gaussian_hmc_reflections_test -tc=cube)
+add_test(NAME ghmc_reflections_test_cross COMMAND gaussian_hmc_reflections_test -tc=cross)
+add_test(NAME ghmc_reflections_test_simplex COMMAND gaussian_hmc_reflections_test -tc=simplex)
+
+
 set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
 
 #set_target_properties(benchmarks_crhmc
@@ -388,3 +394,4 @@ TARGET_LINK_LIBRARIES(logconcave_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT}
 TARGET_LINK_LIBRARIES(crhmc_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT} ${PTHREAD} ${GMP} ${MPSOLVE} ${FFTW3} ${MKL_LINK} QD_LIB coverage_config)
 TARGET_LINK_LIBRARIES(order_polytope lp_solve coverage_config)
 TARGET_LINK_LIBRARIES(matrix_sampling_test lp_solve ${MKL_LINK} coverage_config)
+TARGET_LINK_LIBRARIES(gaussian_hmc_reflections_test lp_solve coverage_config)

test/gaussian_hmc_reflections_test.cpp

@@ -0,0 +1,54 @@
+#include "Eigen/Eigen"
+#include <iostream>
+#include "cartesian_geom/cartesian_kernel.h"
+#include "convex_bodies/hpolytope.h"
+#include "generators/known_polytope_generators.h"
+#include "random_walks/random_walks.hpp"
+#include "misc/misc.h"
+#include "doctest.h"  
+
+typedef double NT;
+typedef Cartesian<NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef BoostRandomNumberGenerator<boost::mt19937, NT, 3> RNGType;
+typedef HPolytope<Point> Polytope;
+
+void test_gaussian_hmc_reflections(const std::string& polytope_type, RNGType& rng, unsigned int walk_length) {
+    Polytope P;
+    if (polytope_type == "cube") {
+        P = generate_cube<Polytope>(3, false, 10.0);
+    } else if (polytope_type == "cross") {
+        P = generate_cross<Polytope>(3, false);
+    } else if (polytope_type == "simplex") {
+        P = generate_simplex<Polytope>(3, false);
+    } else {
+        std::cout << "Unknown polytope type: " << polytope_type << std::endl;
+        return;
+    }
+
+    Point interior_point = P.ComputeInnerBall().first;
+    double L = 1.0; // Step size
+    unsigned int rho = 100; // Max reflections
+    typename GaussianHamiltonianMonteCarloExactWalk::parameters params(L, true, rho, true);
+    typename GaussianHamiltonianMonteCarloExactWalk::Walk<Polytope, RNGType> walk(P, interior_point, 1.0, rng, params);
+
+    walk.apply(P, interior_point, 1.0, walk_length, rng);
+
+    CHECK(P.is_in(interior_point) == -1);  // Use CHECK from doctest
+    std::cout << "Test passed for " << polytope_type << ": Point remains inside the polytope after reflections.\n";
+}
+
+TEST_CASE("GHMC Reflections Test - Cube") {
+    RNGType rng(123456);
+    test_gaussian_hmc_reflections("cube", rng, 10);
+}
+
+TEST_CASE("GHMC Reflections Test - Cross Polytope") {
+    RNGType rng(123456);
+    test_gaussian_hmc_reflections("cross", rng, 10);
+}
+
+TEST_CASE("GHMC Reflections Test - Simplex") {
+    RNGType rng(123456);
+    test_gaussian_hmc_reflections("simplex", rng, 10);
+}