From f06609156716d362a33668e6a3cc2022583f67d2 Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Wed, 7 Aug 2024 13:17:01 +0200
Subject: [PATCH] fix bug with reflection

---
 include/convex_bodies/hpolytope.h              | 18 +++++++-----------
 .../uniform_accelerated_billiard_walk.hpp      | 10 +++++-----
 2 files changed, 12 insertions(+), 16 deletions(-)

diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h
index ccfa9ae1e..f27e7c934 100644
--- a/include/convex_bodies/hpolytope.h
+++ b/include/convex_bodies/hpolytope.h
@@ -1015,17 +1015,13 @@ class HPolytope {
     }
 
     template <typename update_parameters>
-    void compute_reflection_abw(Point &v, Point &p, update_parameters const& params) const {
-            if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) { // MT must be in RowMajor format
-                NT* v_data = v.pointerToData();
-                NT* p_data = p.pointerToData();
-                for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
-                    *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
-                    *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
-                }
-            } else {
-                Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
-                v += a;
+    auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
+         -> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
+            NT* v_data = v.pointerToData();
+            NT* p_data = p.pointerToData();
+            for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
+                *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
+                *(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
             }
     }
 
diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp
index c593b9e6f..604ff202e 100644
--- a/include/random_walks/uniform_accelerated_billiard_walk.hpp
+++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp
@@ -61,7 +61,7 @@ struct AcceleratedBilliardWalk
     {
         typedef typename Polytope::PointType Point;
         typedef typename Polytope::MT MT;
-        typedef typename Polytope::DenseMT DenseMT;
+        typedef typename Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
         typedef typename Point::FT NT;
         using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >; 
         // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
@@ -152,7 +152,7 @@ struct AcceleratedBilliardWalk
                     _p += (_lambda_prev * _v);
                 }
                 T -= _lambda_prev;
-                P.compute_reflection_abw(_v, _p, _update_parameters);
+                P.compute_reflection(_v, _p, _update_parameters);
                 it++;
 
                 while (it < _rho)
@@ -177,7 +177,7 @@ struct AcceleratedBilliardWalk
                         _p += (_lambda_prev * _v);
                     }
                     T -= _lambda_prev;
-                    P.compute_reflection_abw(_v, _p, _update_parameters);
+                    P.compute_reflection(_v, _p, _update_parameters);
                     it++;
                 }
                 _p += _update_parameters.moved_dist * _v;
@@ -298,7 +298,7 @@ struct AcceleratedBilliardWalk
             _lambda_prev = dl * pbpair.first;
             _p += (_lambda_prev * _v);
             T -= _lambda_prev;
-            P.compute_reflection_abw(_v, _p, _update_parameters);
+            P.compute_reflection(_v, _p, _update_parameters);
 
             while (it <= _rho)
             {
@@ -316,7 +316,7 @@ struct AcceleratedBilliardWalk
                 _lambda_prev = dl * pbpair.first;
                 _p += (_lambda_prev * _v);
                 T -= _lambda_prev;
-                P.compute_reflection_abw(_v, _p, _update_parameters);
+                P.compute_reflection(_v, _p, _update_parameters);
                 it++;
             }
         }