From ce1739af2aa00822a320ceb3f5bd567835a02dd8 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Tue, 3 Sep 2024 18:11:26 +0200 Subject: [PATCH 1/3] simplify gabw sampling --- include/convex_bodies/hpolytope.h | 5 ++-- .../gaussian_accelerated_billiard_walk.hpp | 27 ++++++++++++++----- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 9cabcd9ab..f2648aceb 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -1022,14 +1022,15 @@ class HPolytope { } template - NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const ¶ms) const { + NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT &vEv, update_parameters const ¶ms) const { 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 coeff = std::sqrt(vEv / new_vEv); - v = v * coeff; + // v = v * coeff; + vEv = new_vEv; return coeff; } diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index 532d7672e..5176e3a78 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -179,16 +179,20 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; + + T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + T = T / coef; + it++; while (it < _rho) { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - _Av *= coef; - _update_parameters.inner_vi_ak *= coef; - pbpair.first /= coef; + //_Av *= coef; + //_update_parameters.inner_vi_ak *= coef; + //pbpair.first /= coef; if (T <= pbpair.first) { _p += (T * _v); @@ -198,7 +202,11 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; + + T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + T = T / coef; + it++; } if (it == _rho){ @@ -251,15 +259,18 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; + + T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + T = T / coef; while (it <= _rho) { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - _Av *= coef; - _update_parameters.inner_vi_ak *= coef; - pbpair.first /= coef; + //_Av *= coef; + //_update_parameters.inner_vi_ak *= coef; + //pbpair.first /= coef; if (T <= pbpair.first) { _p += (T * _v); @@ -273,7 +284,11 @@ struct GaussianAcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; _p += (_lambda_prev * _v); T -= _lambda_prev; + + T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); + T = T / coef; + it++; } } From 19565dc1eacb2dd43b2f69f042675099b5549710 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 6 Sep 2024 14:05:57 +0200 Subject: [PATCH 2/3] optimize gabw for sparse polytopes --- include/convex_bodies/hpolytope.h | 40 +++++-- .../gaussian_accelerated_billiard_walk.hpp | 102 ++++++++++++------ .../uniform_accelerated_billiard_walk.hpp | 2 +- test/sampling_test.cpp | 22 +++- 4 files changed, 119 insertions(+), 47 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index f2648aceb..03be9bbbc 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -1010,9 +1010,10 @@ class HPolytope { // 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 template auto compute_reflection(Point &v, Point &p, update_parameters const& params) const - -> std::enable_if_t> && !std::is_same_v, void> { // MT must be in RowMajor format + -> std::enable_if_t> && !std::is_same_v, void> { NT* v_data = v.pointerToData(); NT* p_data = p.pointerToData(); for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { @@ -1021,15 +1022,38 @@ class HPolytope { } } - template - NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT &vEv, update_parameters const ¶ms) const { + // 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 { + + NT new_vEv; + if constexpr (!std::is_same_v>) { + 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>) { + 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::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::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; } diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index 5176e3a78..b2d1b9644 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -47,12 +47,13 @@ struct GaussianAcceleratedBilliardWalk struct update_parameters { update_parameters() - : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0) + : facet_prev(0), hit_ball(false), inner_vi_ak(0.0), ball_inner_norm(0.0), moved_dist(0.0) {} int facet_prev; bool hit_ball; double inner_vi_ak; double ball_inner_norm; + double moved_dist; }; parameters param; @@ -67,9 +68,15 @@ struct GaussianAcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Polytope::DenseMT DenseMT; + typedef typename Polytope::MT MT; + 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 >; + // ( ( ( )) ( ( ) ) ( ) ) + // AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise void computeLcov(const E_type E) { @@ -110,7 +117,11 @@ struct GaussianAcceleratedBilliardWalk _L = compute_diameter::template compute(P); computeLcov(E); _E = E; - _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + if constexpr (std::is_same>::value) { + _AA = (P.get_mat() * P.get_mat().transpose()); + } else { + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -131,7 +142,11 @@ struct GaussianAcceleratedBilliardWalk ::template compute(P); computeLcov(E); _E = E; - _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + if constexpr (std::is_same>::value) { + _AA = (P.get_mat() * P.get_mat().transpose()); + } else { + _AA.noalias() = (DenseMT)(P.get_mat() * P.get_mat().transpose()); + } _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental) initialize(P, p, rng); } @@ -148,6 +163,12 @@ 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 pbpair; - if(!was_reset) { - pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _update_parameters); - } else { - pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); - was_reset = false; - } + std::pair pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters); if (T <= pbpair.first) { _p += (T * _v); @@ -177,7 +192,18 @@ struct GaussianAcceleratedBilliardWalk } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist = _lambda_prev; + NT* Ar_data = _lambdas.data(); + NT* Av_data = _Av.data(); + for(int i = 0; i < P.num_of_hyperplanes(); ++i) { + _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + } + // rebuild the heap with the new values of (b - Ar) / Av + _distances_set.rebuild(_update_parameters.moved_dist); + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; T = T * coef; @@ -188,11 +214,14 @@ struct GaussianAcceleratedBilliardWalk while (it < _rho) { - std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - //_Av *= coef; - //_update_parameters.inner_vi_ak *= coef; - //pbpair.first /= coef; + std::pair pbpair; + if constexpr (std::is_same>::value) { + pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, + _distances_set, _AA, _update_parameters); + } else { + pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, + _AA, _update_parameters); + } if (T <= pbpair.first) { _p += (T * _v); @@ -200,7 +229,11 @@ struct GaussianAcceleratedBilliardWalk break; } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist += _lambda_prev; + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; T = T * coef; @@ -209,9 +242,10 @@ struct GaussianAcceleratedBilliardWalk it++; } + _p += _update_parameters.moved_dist * _v; + _update_parameters.moved_dist = 0.0; if (it == _rho){ _p = p0; - was_reset = true; } } @@ -227,21 +261,24 @@ struct GaussianAcceleratedBilliardWalk { unsigned int n = P.dimension(); const NT dl = 0.995; - was_reset = false; _lambdas.setZero(P.num_of_hyperplanes()); _Av.setZero(P.num_of_hyperplanes()); _p = p; - _AE.noalias() = (DenseMT)(P.get_mat() * _E); - _AEA = _AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum(); - /* - _AEA.resize(P.num_of_hyperplanes()); - for(int i = 0; i < P.num_of_hyperplanes(); ++i) { - _AEA(i) = _AE.row(i).dot(P.get_mat().row(i)); - }*/ + DenseMT temp_AE; + temp_AE.noalias() = (DenseMT)(P.get_mat() * _E); + _AEA = temp_AE.cwiseProduct((DenseMT)P.get_mat()).rowwise().sum(); + if constexpr (std::is_same_v) + _AE = temp_AE; + else + _AE = temp_AE.sparseView(); + } + _distances_set = BoundaryOracleHeap(P.num_of_hyperplanes()); + _v = GetDirection::apply(n, rng, false); _v = Point(_L_cov.template triangularView() * _v.getCoefficients()); + NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; @@ -268,9 +305,6 @@ struct GaussianAcceleratedBilliardWalk { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); - //_Av *= coef; - //_update_parameters.inner_vi_ak *= coef; - //pbpair.first /= coef; if (T <= pbpair.first) { _p += (T * _v); @@ -288,7 +322,7 @@ struct GaussianAcceleratedBilliardWalk T = T * coef; coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters); T = T / coef; - + it++; } } @@ -297,16 +331,16 @@ struct GaussianAcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - DenseMT _AA; + AA_type _AA; E_type _L_cov; // LL' = inv(E) - DenseMT _AE; + AE_type _AE; E_type _E; VT _AEA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - bool was_reset; + BoundaryOracleHeap _distances_set; }; }; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index f91abee4b..77660c745 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -275,7 +275,7 @@ struct AcceleratedBilliardWalk it++; while (it < _rho) - { + { std::pair pbpair; if constexpr (std::is_same>::value) { pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index a15f40102..ba71d1398 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -320,8 +320,11 @@ void call_test_gabw(){ Point StartingPoint(d); std::list randPoints; + std::chrono::time_point start, stop; + start = std::chrono::high_resolution_clock::now(); + std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Skinny-H-cube10" << std::endl; - P = generate_skinny_cube(10); + P = generate_skinny_cube(d); Point p = P.ComputeInnerBall().first; @@ -340,6 +343,11 @@ void call_test_gabw(){ RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints, push_back_policy, rng); + stop = std::chrono::high_resolution_clock::now(); + + std::chrono::duration total_time = stop - start; + std::cout << "Done in " << total_time.count() << '\n'; + MT samples(d, numpoints); unsigned int jj = 0; for (typename std::list::iterator rpit = randPoints.begin(); rpit != randPoints.end(); rpit++, jj++) @@ -352,11 +360,11 @@ void call_test_gabw(){ RNGType Srng(d); typedef Eigen::SparseMatrix SparseMT; - typedef HPolytope SparseHpoly; + typedef HPolytope> SparseHpoly; std::list Points; SparseHpoly SP; - SP = generate_skinny_cube(10); + SP = generate_skinny_cube(d); std::cout << "--- Testing Gaussian Accelerated Billiard Walk for Sparse Skinny-H-cube10" << std::endl; @@ -371,15 +379,21 @@ void call_test_gabw(){ > sparsewalk; typedef MultivariateGaussianRandomPointGenerator SparseRandomPointGenerator; + start = std::chrono::high_resolution_clock::now(); ellipsoid = compute_inscribed_ellipsoid - (SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); + ((SparseMT)SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0)); const SparseMT SE = get<0>(ellipsoid).sparseView(); SparseRandomPointGenerator::apply(SP, p, SE, numpoints, 1, Points, push_back_policy, Srng); + stop = std::chrono::high_resolution_clock::now(); + + total_time = stop - start; + std::cout << "Done in " << total_time.count() << '\n'; + jj = 0; MT Ssamples(d, numpoints); for (typename std::list::iterator rpit = Points.begin(); rpit != Points.end(); rpit++, jj++) From 4fc2b3fd7e6b8dd7f1bb54c468ad3f4d7c74f6d9 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Sat, 14 Sep 2024 17:50:36 +0200 Subject: [PATCH 3/3] minor changes --- include/convex_bodies/hpolytope.h | 6 +- .../accelerated_billiard_walk_utils.hpp | 133 ++++++++++++++++++ .../gaussian_accelerated_billiard_walk.hpp | 24 ++-- .../uniform_accelerated_billiard_walk.hpp | 120 +--------------- 4 files changed, 148 insertions(+), 135 deletions(-) create mode 100644 include/random_walks/accelerated_billiard_walk_utils.hpp 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