fix bug with reflection

GeomScale · Aug 7, 2024 · f066091 · f066091
1 parent c61c4bd
commit f066091
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 16 deletions.
diff --git 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
@@ -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++;
             }
         }