compute_reflection_abw_sparse and improve comments

Author: lucaperju

.gitignore

@@ -9,4 +9,5 @@ test/Testing/Temporary/CTestCostData.txt
 build/
 .vscode
 .DS_Store
-.cache
+.cache
+shell.nix

include/convex_bodies/hpolytope.h

@@ -580,17 +580,23 @@ class HPolytope {
         NT inner_prev = params.inner_vi_ak;
         NT* Av_data = Av.data();
 
+        // Updating Av due to the change in direction caused by the previous reflection
         // Av += (-2.0 * inner_prev) * AA.col(params.facet_prev)
 
         for (Eigen::SparseMatrix<double>::InnerIterator it(AA, params.facet_prev); it; ++it) {
 
+            // val(row) - params.moved_dist = The distance until we would hit the facet given by row
+            // all those values are stored inside distances_set for quick retrieval of the minimum
 
+            // Before updating Av(row)
             // val(row) = (b(row) - Ar(row)) / Av(row) + params.moved_dist
             // (val(row) - params.moved_dist) = (b(row) - Ar(row)) / Av(row)
             // (val(row) - params.moved_dist) * Av(row) = b(row) - Ar(row)
+            // b(row) - Ar(row) = (val(row) - params.moved_dist) * Av(row)
 
             *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value();
 
+            // After updating Av(row)
             // b(row) - Ar(row) = (old_val(row) - params.moved_dist) * old_Av(row) 
             // new_val(row) = (b(row) - Ar(row)                                ) / new_Av(row) + params.moved_dist 
             // new_val(row) = ((old_val(row) - params.moved_dist) * old_Av(row)) / new_Av(row) + params.moved_dist
@@ -608,8 +614,12 @@ class HPolytope {
             distances_set.change_val(it.row(), val, params.moved_dist);
         }
 
+        // Finding the distance to the closest facet and its row
         std::pair<NT, int> ans = distances_set.get_min();
+
+        // Subtract params.moved_dist to obtain the actual distance to the facet
         ans.first -= params.moved_dist;
+
         params.inner_vi_ak = *(Av_data + ans.second);
         params.facet_prev = ans.second;
 
@@ -1002,24 +1012,24 @@ class HPolytope {
       return total;
     }
 
+    // Updates the velocity vector v and the position vector p after a reflection
     template <typename update_parameters>
     void compute_reflection(Point &v, Point const&, update_parameters const& params) const {
             Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
             v += a;
     }
 
-    // updates the velocity vector v and the position vector p after a reflection
-    // the real value of p is given by p + moved_dist * v
-    // MT must be sparse, in RowMajor format
+    // Only to be called when MT is in RowMajor format
+    // The real value of p is given by p + params.moved_dist * v
     template <typename update_parameters>
-    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> {
-            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();
-            }
+    auto compute_reflection_abw_sparse(Point &v, Point &p, update_parameters const& params) const
+    {
+        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();
+        }
     }
 
     // function to compute reflection for GaBW random walk

include/random_walks/uniform_accelerated_billiard_walk.hpp

@@ -155,7 +155,11 @@ struct AcceleratedBilliardWalk
                     _p += (_lambda_prev * _v);
                 }
                 T -= _lambda_prev;
-                P.compute_reflection_abw(_v, _p, _update_parameters);
+                if constexpr (SPARSE) {
+                    P.compute_reflection_abw_sparse(_v, _p, _update_parameters);
+                } else {
+                    P.compute_reflection(_v, _p, _update_parameters);
+                }
                 it++;
 
                 while (it < _rho)
@@ -180,7 +184,11 @@ struct AcceleratedBilliardWalk
                         _p += (_lambda_prev * _v);
                     }
                     T -= _lambda_prev;
-                    P.compute_reflection_abw(_v, _p, _update_parameters);
+                    if constexpr (SPARSE) {
+                        P.compute_reflection_abw_sparse(_v, _p, _update_parameters);
+                    } else {
+                        P.compute_reflection(_v, _p, _update_parameters);
+                    }
                     it++;
                 }
                 _p += _update_parameters.moved_dist * _v;
@@ -302,7 +310,11 @@ struct AcceleratedBilliardWalk
             _lambda_prev = dl * pbpair.first;
             _p += (_lambda_prev * _v);
             T -= _lambda_prev;
-            P.compute_reflection_abw(_v, _p, _update_parameters);
+            if constexpr (SPARSE) {
+                P.compute_reflection_abw_sparse(_v, _p, _update_parameters);
+            } else {
+                P.compute_reflection(_v, _p, _update_parameters);
+            }
 
             while (it <= _rho)
             {
@@ -320,7 +332,11 @@ struct AcceleratedBilliardWalk
                 _lambda_prev = dl * pbpair.first;
                 _p += (_lambda_prev * _v);
                 T -= _lambda_prev;
-                P.compute_reflection_abw(_v, _p, _update_parameters);
+                if constexpr (SPARSE) {
+                    P.compute_reflection_abw_sparse(_v, _p, _update_parameters);
+                } else {
+                    P.compute_reflection(_v, _p, _update_parameters);
+                }
                 it++;
             }
         }