diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 03be9bbbc..24cfd7284 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -1024,8 +1024,8 @@ class HPolytope { // function to compute reflection for GaBW random walk // compatible when the polytope is both dense or sparse - template - NT compute_reflection(Point &v, Point &p, AE_type const &AE, VT const &AEA, NT &vEv, update_parameters const ¶ms) const { + template + NT compute_reflection(Point &v, Point &p, NT &vEv, DenseSparseMT const &AE, VT const &AEA, update_parameters const ¶ms) const { NT new_vEv; if constexpr (!std::is_same_v>) { @@ -1035,7 +1035,7 @@ class HPolytope { v += a; } else { - if constexpr(!std::is_same_v>) { + if constexpr(!std::is_same_v>) { 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 { diff --git a/include/random_walks/accelerated_billiard_walk_utils.hpp b/include/random_walks/accelerated_billiard_walk_utils.hpp new file mode 100644 index 000000000..24cf32d76 --- /dev/null +++ b/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 +#include + +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 +class BoundaryOracleHeap { +public: + int n, heap_size; + std::vector> heap; + std::vector> 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 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 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 diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index b2d1b9644..a8a7d334e 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -21,6 +21,7 @@ #include "generators/boost_random_number_generator.hpp" #include "sampling/ellipsoid.hpp" #include "random_walks/uniform_billiard_walk.hpp" +#include "random_walks/accelerated_billiard_walk_utils.hpp" #include "random_walks/compute_diameter.hpp" @@ -72,11 +73,10 @@ struct GaussianAcceleratedBilliardWalk typedef typename Eigen::Matrix DenseMT; typedef typename Polytope::VT VT; typedef typename Point::FT NT; - using AA_type = std::conditional_t< std::is_same_v>, typename Eigen::SparseMatrix, DenseMT >; // AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise - using AE_type = std::conditional_t< std::is_same_v> && std::is_base_of, E_type>::value, typename Eigen::SparseMatrix, DenseMT >; - // ( ( ( )) ( ( ) ) ( ) ) + using AA_type = std::conditional_t< std::is_same_v>, typename Eigen::SparseMatrix, DenseMT >; // AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise + using AE_type = std::conditional_t< std::is_same_v> && std::is_base_of, E_type>::value, typename Eigen::SparseMatrix, DenseMT >; void computeLcov(const E_type E) { @@ -163,12 +163,6 @@ struct GaussianAcceleratedBilliardWalk int it; NT coef = 1.0; NT vEv; - typename Point::Coeff b; - NT* b_data; - if constexpr (std::is_same>::value) { - b = P.get_vec(); - b_data = b.data(); - } for (auto j=0u; j>::value) { + typename Point::Coeff b; + NT* b_data; + b = P.get_vec(); + b_data = b.data(); _update_parameters.moved_dist = _lambda_prev; NT* Ar_data = _lambdas.data(); NT* Av_data = _Av.data(); @@ -207,7 +205,7 @@ struct GaussianAcceleratedBilliardWalk T -= _lambda_prev; T = T * coef; - coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters); T = T / coef; it++; @@ -237,7 +235,7 @@ struct GaussianAcceleratedBilliardWalk T -= _lambda_prev; T = T * coef; - coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters); T = T / coef; it++; @@ -298,7 +296,7 @@ struct GaussianAcceleratedBilliardWalk T -= _lambda_prev; T = T * coef; - coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters); T = T / coef; while (it <= _rho) @@ -320,7 +318,7 @@ struct GaussianAcceleratedBilliardWalk T -= _lambda_prev; T = T * coef; - coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters); T = T / coef; it++; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 77660c745..f68a21a2c 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -13,125 +13,7 @@ #include "sampling/sphere.hpp" #include -#include -#include - -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 -class BoundaryOracleHeap { -public: - int n, heap_size; - std::vector> heap; - std::vector> 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 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 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); - } - } - } -}; +#include "random_walks/accelerated_billiard_walk_utils.hpp" // Billiard walk which accelarates each step for uniform distribution