Sparse optimizations for GaBW sampling

.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,34 +1012,59 @@ 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
+    // 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> { // 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();
-            }
+    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();
+        }
     }
 
-    template <typename update_parameters>
-    NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) const {
+    // function to compute reflection for GaBW random walk
+    // compatible when the polytope is both dense or sparse
+    template <typename DenseSparseMT, typename update_parameters>
+    NT compute_reflection(Point &v, Point &p, NT &vEv, DenseSparseMT const &AE, VT const &AEA, update_parameters const &params) const {
+
+            NT new_vEv;
+            if constexpr (!std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
+                Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
+                VT x = v.getCoefficients();
+                new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
+                v += a;
+            } else {
+
+                if constexpr(!std::is_same_v<DenseSparseMT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
+                    VT x = v.getCoefficients();
+                    new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
+                } else {
+                    new_vEv = vEv + 4.0 * params.inner_vi_ak * params.inner_vi_ak * AEA(params.facet_prev);
+                    NT* v_data = v.pointerToData();
+                    for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(AE, params.facet_prev); it; ++it) {
+                        new_vEv -= 4.0 * params.inner_vi_ak * it.value() * *(v_data + it.col());
+                    }
+                }
 
-            Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
-            VT x = v.getCoefficients();
-            NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
-            v += a;
+                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();
+                }
+            }
             NT coeff = std::sqrt(vEv / new_vEv);
-            v = v * coeff;
+            vEv = new_vEv;
             return coeff;
     }
 

include/random_walks/accelerated_billiard_walk_utils.hpp

@@ -0,0 +1,133 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2024 Luca Perju
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef ACCELERATED_BILLIARD_WALK_UTILS_HPP
+#define ACCELERATED_BILLIARD_WALK_UTILS_HPP
+
+#include <Eigen/Eigen>
+#include <vector>
+
+const double eps = 1e-10;
+
+// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it
+// vec[i].first contains the value of (b(i) - Ar(i))/Av(i) + moved_dist, where moved_dist is the total distance that the point has travelled so far
+// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive)
+template<typename NT>
+class BoundaryOracleHeap {
+public:
+    int n, heap_size;
+    std::vector<std::pair<NT, int>> heap;
+    std::vector<std::pair<NT, int>> vec;
+
+private:
+    int siftDown(int index) {
+        while((index << 1) + 1 < heap_size) {
+            int child = (index << 1) + 1;
+            if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) {
+                child += 1;
+            }
+            if(heap[child].first < heap[index].first - eps)
+            {
+                std::swap(heap[child], heap[index]);
+                std::swap(vec[heap[child].second].second, vec[heap[index].second].second);
+                index = child;
+            } else {
+                return index;
+            }
+        }
+        return index;
+    }
+
+    int siftUp(int index) {
+        while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) {
+            std::swap(heap[(index - 1) >> 1], heap[index]);
+            std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second);
+            index = (index - 1) >> 1;
+        }
+        return index;
+    }
+
+    // takes the index of a facet, and (in case it is in the heap) removes it from the heap.
+    void remove (int index) {
+        index = vec[index].second;
+        if(index == -1) {
+            return;
+        }
+        std::swap(heap[heap_size - 1], heap[index]);
+        std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second);
+        vec[heap[heap_size - 1].second].second = -1;
+        heap_size -= 1;
+        index = siftDown(index);
+        siftUp(index);
+    }
+
+    // inserts a new value into the heap, with its associated facet
+    void insert (const std::pair<NT, int> val) {
+        vec[val.second].second = heap_size;
+        vec[val.second].first = val.first;
+        heap[heap_size++] = val;
+        siftUp(heap_size - 1);
+    }
+
+public:
+    BoundaryOracleHeap() {}
+
+    BoundaryOracleHeap(int n) : n(n), heap_size(0) {
+        heap.resize(n);
+        vec.resize(n);
+    }
+
+    // rebuilds the heap with the existing values from vec
+    // O(n)
+    void rebuild (const NT &moved_dist) {
+        heap_size = 0;
+        for(int i = 0; i < n; ++i) {
+            vec[i].second = -1;
+            if(vec[i].first - eps > moved_dist) {
+                vec[i].second = heap_size;
+                heap[heap_size++] = {vec[i].first, i};
+            }
+        }
+        for(int i = heap_size - 1; i >= 0; --i) {
+            siftDown(i);
+        }
+    }
+
+    // returns (b(i) - Ar(i))/Av(i) + moved_dist
+    // O(1)
+    NT get_val (const int &index) {
+        return vec[index].first;
+    }
+
+    // returns the nearest facet
+    // O(1)
+    std::pair<NT, int> get_min () {
+        return heap[0];
+    }
+
+    // changes the stored value for a given facet, and updates the heap accordingly
+    // O(logn)
+    void change_val(const int& index, const NT& new_val, const NT& moved_dist) {
+        if(new_val < moved_dist - eps) {
+            vec[index].first = new_val;
+            remove(index);
+        } else {
+            if(vec[index].second == -1) {
+                insert({new_val, index});
+            } else {
+                heap[vec[index].second].first = new_val;
+                vec[index].first = new_val;
+                siftDown(vec[index].second);
+                siftUp(vec[index].second);
+            }
+        }
+    }
+};
+
+
+#endif