From ae5b45f9dcff63b2e958571a66625e41ee1a8019 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 5 Aug 2024 13:59:31 +0200 Subject: [PATCH 01/15] support for eigen sparse rowmajor matrix --- include/convex_bodies/hpolytope.h | 14 +++++++++---- .../preprocess/rounding_util_functions.hpp | 20 +++++++++---------- .../gaussian_accelerated_billiard_walk.hpp | 2 +- .../uniform_accelerated_billiard_walk.hpp | 8 ++++---- test/sampling_test.cpp | 2 +- 5 files changed, 26 insertions(+), 20 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 5f37be325..ce5d810f3 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -525,7 +525,7 @@ class HPolytope { VT& Av, NT const& lambda_prev, DenseMT const& AA, - update_parameters& params) const + update_parameters& params) const // { NT min_plus = std::numeric_limits::max(); @@ -954,9 +954,15 @@ class HPolytope { template void compute_reflection(Point &v, const Point &, update_parameters const& params) const { - - Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); - v += a; + if constexpr (std::is_same>::value) { // is faster only if MT is RowMajor + NT* v_data = v.pointerToData(); + for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { + *(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value(); + } + } else { + Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); + v += a; + } } template diff --git a/include/preprocess/rounding_util_functions.hpp b/include/preprocess/rounding_util_functions.hpp index a189890e3..bb07d5b30 100644 --- a/include/preprocess/rounding_util_functions.hpp +++ b/include/preprocess/rounding_util_functions.hpp @@ -57,7 +57,7 @@ initialize_chol(MT const& mat) if constexpr (std::is_same::value) { return std::make_unique>(); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { auto llt = std::make_unique>(); llt->analyzePattern(mat); @@ -78,7 +78,7 @@ initialize_chol(MT const& A_trans, MT const& A) if constexpr (std::is_same::value) { return std::make_unique>(); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { MT mat = A_trans * A; return initialize_chol(mat); @@ -99,7 +99,7 @@ inline static VT solve_vec(std::unique_ptr const& llt, { llt->compute(H); return llt->solve(b); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { llt->factorize(H); return llt->solve(b); @@ -122,7 +122,7 @@ solve_mat(std::unique_ptr const& llt, llt->compute(H); logdetE = llt->matrixL().toDenseMatrix().diagonal().array().log().sum(); return llt->solve(mat); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { llt->factorize(H); logdetE = llt->matrixL().nestedExpression().diagonal().array().log().sum(); @@ -143,7 +143,7 @@ inline static void update_Atrans_Diag_A(MT &H, MT const& A_trans, if constexpr (std::is_same::value) { H.noalias() = A_trans * D * A; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = A_trans * D * A; } else @@ -161,7 +161,7 @@ inline static void update_Diag_A(MT &H, diag_MT const& D, MT const& A) if constexpr (std::is_same::value) { H.noalias() = D * A; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = D * A; } else @@ -179,7 +179,7 @@ inline static void update_A_Diag(MT &H, MT const& A, diag_MT const& D) if constexpr (std::is_same::value) { H.noalias() = A * D; - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { H = A * D; } else @@ -198,7 +198,7 @@ get_mat_prod_op(MT const& E) if constexpr (std::is_same::value) { return std::make_unique>(E); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { return std::make_unique>(E); } else @@ -249,7 +249,7 @@ init_Bmat(MT &B, int const n, MT const& A_trans, MT const& A) if constexpr (std::is_same::value) { B.resize(n+1, n+1); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { // Initialize the structure of matrix B typedef Eigen::Triplet triplet; @@ -299,7 +299,7 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d, B.block(n, 0, 1, n).noalias() = AtDe.transpose(); B(n, n) = d.sum(); B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1); - } else if constexpr (std::is_same::value) + } else if constexpr (std::is_base_of, MT >::value) { MT AtD_A = AtD * A; int k = 0; diff --git a/include/random_walks/gaussian_accelerated_billiard_walk.hpp b/include/random_walks/gaussian_accelerated_billiard_walk.hpp index d4cd7d9df..532d7672e 100644 --- a/include/random_walks/gaussian_accelerated_billiard_walk.hpp +++ b/include/random_walks/gaussian_accelerated_billiard_walk.hpp @@ -73,7 +73,7 @@ struct GaussianAcceleratedBilliardWalk void computeLcov(const E_type E) { - if constexpr (std::is_same >::value) { + if constexpr (std::is_base_of, E_type >::value) { Eigen::SimplicialLLT lltofE; lltofE.compute(E); if (lltofE.info() != Eigen::Success) { diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 60ab38263..eff36189e 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -136,16 +136,16 @@ struct AcceleratedBilliardWalk { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, - _AA, _update_parameters); + _AA, _update_parameters); // if (T <= pbpair.first) { - _p += (T * _v); + _p += (T * _v);// _lambda_prev = T; break; } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + _p += (_lambda_prev * _v);// T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + P.compute_reflection(_v, _p, _update_parameters);// it++; } if (it == _rho) { diff --git a/test/sampling_test.cpp b/test/sampling_test.cpp index a03e46684..a15f40102 100644 --- a/test/sampling_test.cpp +++ b/test/sampling_test.cpp @@ -399,7 +399,7 @@ template void call_test_sparse(){ typedef Cartesian Kernel; typedef typename Kernel::Point Point; - typedef HPolytope> Hpolytope; + typedef HPolytope> Hpolytope; typedef typename Hpolytope::MT MT; typedef typename Hpolytope::DenseMT DenseMT; typedef Eigen::Matrix VT; From 37f0ca08f459ee9b0c237ee0d1e3562e398190ba Mon Sep 17 00:00:00 2001 From: lucaperju Date: Mon, 5 Aug 2024 18:27:30 +0200 Subject: [PATCH 02/15] optimize updating the position in abw --- include/convex_bodies/hpolytope.h | 10 ++++--- .../uniform_accelerated_billiard_walk.hpp | 26 +++++++++++++------ 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index ce5d810f3..ff7d3347d 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -525,7 +525,7 @@ class HPolytope { VT& Av, NT const& lambda_prev, DenseMT const& AA, - update_parameters& params) const // + update_parameters& params) const// { NT min_plus = std::numeric_limits::max(); @@ -918,7 +918,7 @@ class HPolytope { normalized = true; } - void compute_reflection(Point& v, Point const&, int const& facet) const + void compute_reflection(Point& v, Point& p, int const& facet) const { v += -2 * v.dot(A.row(facet)) * A.row(facet); } @@ -953,11 +953,13 @@ class HPolytope { } template - void compute_reflection(Point &v, const Point &, update_parameters const& params) const { - if constexpr (std::is_same>::value) { // is faster only if MT is RowMajor + void compute_reflection(Point &v, Point &p, update_parameters const& params) const { + if constexpr (std::is_same>::value) { // is faster only if MT is in RowMajor format 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(); } } else { Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev)); diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index eff36189e..aefc9fa7d 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -39,12 +39,13 @@ struct AcceleratedBilliardWalk 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; @@ -58,7 +59,8 @@ struct AcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Eigen::Matrix MT; + typedef typename Polytope::MT MT; + typedef typename Polytope::DenseMT DenseMT; typedef typename Point::FT NT; template @@ -70,7 +72,7 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); - _AA.noalias() = (MT)(P.get_mat() * P.get_mat().transpose()); + _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); } @@ -86,7 +88,7 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - _AA.noalias() = (MT)(P.get_mat() * P.get_mat().transpose()); + _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); } @@ -132,22 +134,30 @@ struct AcceleratedBilliardWalk P.compute_reflection(_v, _p, _update_parameters); it++; + _update_parameters.moved_dist = 0.0; + while (it < _rho) { std::pair pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, _AA, _update_parameters); // if (T <= pbpair.first) { - _p += (T * _v);// + _p += (T * _v); _lambda_prev = T; 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);// done + } T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters);// + P.compute_reflection(_v, _p, _update_parameters);// done it++; } + _p += _update_parameters.moved_dist * _v; + _update_parameters.moved_dist = 0.0; if (it == _rho) { _p = p0; was_reset = true; @@ -312,7 +322,7 @@ struct AcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - MT _AA; + DenseMT _AA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; From 725606b2c6d5008928e46b9ecd1a6ef0384e88e1 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 7 Aug 2024 11:09:26 +0200 Subject: [PATCH 03/15] finish sparse optimizations for UniformABW --- include/convex_bodies/hpolytope.h | 60 +++++++++++++++++- .../uniform_accelerated_billiard_walk.hpp | 63 +++++++++++++------ 2 files changed, 101 insertions(+), 22 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index ff7d3347d..6f9e7eaba 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -525,7 +525,7 @@ class HPolytope { VT& Av, NT const& lambda_prev, DenseMT const& AA, - update_parameters& params) const// + update_parameters& params) const { NT min_plus = std::numeric_limits::max(); @@ -567,6 +567,62 @@ class HPolytope { } + + template + std::pair line_positive_intersect(Point const& r, + Point const& v, + VT& Ar, + VT& Av, + VT& distances_vec, + set_type& distances_set, + AA_type const& AA, + update_parameters& params) const + { + NT inner_prev = params.inner_vi_ak; + + // real Ar = Ar + params.moved_dist * Av + // real r = r + params.moved_dist * v + // real distances_vec = distances_vec - params.moved_dist + + NT* Av_data = Av.data(); + NT* Ar_data = Ar.data(); + NT* dvec_data = distances_vec.data(); + const NT* b_data = b.data(); + for (Eigen::SparseMatrix::InnerIterator it(AA, params.facet_prev); it; ++it) { + + *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); + *(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value(); + + distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); + + *(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row()); + + if(*(dvec_data + it.row()) > params.moved_dist) + distances_set.insert(std::make_pair(*(dvec_data + it.row()), it.row())); + } + + auto it = distances_set.upper_bound(std::make_pair(params.moved_dist, 0)); + + if(it == distances_set.end()) { + std::cout << "something went wrong when trying to get lowest positive value" << std::endl; + throw "all values from the set were negative"; + } + + std::pair ans = (*it); + ans.first -= params.moved_dist; + + params.inner_vi_ak = *(Av_data + ans.second); + params.facet_prev = ans.second; + + /*if(ans.first < 0.00000001) { + std::cout << "distance of 0 found" << std::endl; + exit(0); + }*/ + + return ans; + } + + template std::pair line_positive_intersect(Point const& r, Point const& v, @@ -954,7 +1010,7 @@ class HPolytope { template void compute_reflection(Point &v, Point &p, update_parameters const& params) const { - if constexpr (std::is_same>::value) { // is faster only if MT is in RowMajor format + if constexpr (std::is_same>::value) { // MT must be in RowMajor format NT* v_data = v.pointerToData(); NT* p_data = p.pointerToData(); for(Eigen::SparseMatrix::InnerIterator it(A, params.facet_prev); it; ++it) { diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index aefc9fa7d..b0ee860a7 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -13,6 +13,7 @@ #include "sampling/sphere.hpp" #include +#include // Billiard walk which accelarates each step for uniform distribution @@ -62,6 +63,8 @@ struct AcceleratedBilliardWalk typedef typename Polytope::MT MT; typedef typename Polytope::DenseMT DenseMT; 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 template Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng) @@ -72,7 +75,11 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); + 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); } @@ -88,7 +95,11 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); + 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); } @@ -114,13 +125,7 @@ struct AcceleratedBilliardWalk Point p0 = _p; it = 0; - std::pair 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); @@ -129,18 +134,37 @@ struct AcceleratedBilliardWalk } _lambda_prev = dl * pbpair.first; - _p += (_lambda_prev * _v); + if constexpr (std::is_same>::value) { + _update_parameters.moved_dist = _lambda_prev; + _distances_set.clear(); + _distances_vec.setZero(P.num_of_hyperplanes()); + typename Point::Coeff b = P.get_vec(); + NT* b_data = b.data(); + NT* dvec_data = _distances_vec.data(); + NT* Ar_data = _lambdas.data(); + NT* Av_data = _Av.data(); + for(int i = 0; i < P.num_of_hyperplanes(); ++i) { + *(dvec_data + i) = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + if(*(dvec_data + i) > _update_parameters.moved_dist) + _distances_set.insert(std::make_pair(*(dvec_data + i), i)); + } + } else { + _p += (_lambda_prev * _v); + } T -= _lambda_prev; P.compute_reflection(_v, _p, _update_parameters); it++; - _update_parameters.moved_dist = 0.0; - while (it < _rho) - { - std::pair pbpair - = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, - _AA, _update_parameters); // + { + std::pair pbpair; + if constexpr (std::is_same>::value) { + pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _distances_vec, + _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); _lambda_prev = T; @@ -150,17 +174,16 @@ struct AcceleratedBilliardWalk if constexpr (std::is_same>::value) { _update_parameters.moved_dist += _lambda_prev; } else { - _p += (_lambda_prev * _v);// done + _p += (_lambda_prev * _v); } T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters);// done + P.compute_reflection(_v, _p, _update_parameters); it++; } _p += _update_parameters.moved_dist * _v; _update_parameters.moved_dist = 0.0; if (it == _rho) { _p = p0; - was_reset = true; } } p = _p; @@ -256,7 +279,6 @@ struct AcceleratedBilliardWalk { 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; @@ -322,12 +344,13 @@ struct AcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - DenseMT _AA; + AA_type _AA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - bool was_reset; + typename Point::Coeff _distances_vec; + std::set> _distances_set; }; }; From c61c4bdaa77cc3a228438c83b335f13dac5fbb9b Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 7 Aug 2024 11:37:07 +0200 Subject: [PATCH 04/15] temporary fix to failing tests --- include/convex_bodies/hpolytope.h | 10 ++++++++-- .../random_walks/uniform_accelerated_billiard_walk.hpp | 8 ++++---- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 6f9e7eaba..ccfa9ae1e 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -974,7 +974,7 @@ class HPolytope { normalized = true; } - void compute_reflection(Point& v, Point& p, int const& facet) const + void compute_reflection(Point& v, Point const&, int const& facet) const { v += -2 * v.dot(A.row(facet)) * A.row(facet); } @@ -1009,7 +1009,13 @@ class HPolytope { } template - void compute_reflection(Point &v, Point &p, update_parameters const& params) const { + 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; + } + + template + void compute_reflection_abw(Point &v, Point &p, update_parameters const& params) const { if constexpr (std::is_same>::value) { // MT must be in RowMajor format NT* v_data = v.pointerToData(); NT* p_data = p.pointerToData(); diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index b0ee860a7..c593b9e6f 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -152,7 +152,7 @@ struct AcceleratedBilliardWalk _p += (_lambda_prev * _v); } T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + P.compute_reflection_abw(_v, _p, _update_parameters); it++; while (it < _rho) @@ -177,7 +177,7 @@ struct AcceleratedBilliardWalk _p += (_lambda_prev * _v); } T -= _lambda_prev; - P.compute_reflection(_v, _p, _update_parameters); + P.compute_reflection_abw(_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(_v, _p, _update_parameters); + P.compute_reflection_abw(_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(_v, _p, _update_parameters); + P.compute_reflection_abw(_v, _p, _update_parameters); it++; } } From f06609156716d362a33668e6a3cc2022583f67d2 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 7 Aug 2024 13:17:01 +0200 Subject: [PATCH 05/15] fix bug with reflection --- include/convex_bodies/hpolytope.h | 18 +++++++----------- .../uniform_accelerated_billiard_walk.hpp | 10 +++++----- 2 files changed, 12 insertions(+), 16 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index ccfa9ae1e..f27e7c934 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -1015,17 +1015,13 @@ class HPolytope { } template - void compute_reflection_abw(Point &v, Point &p, update_parameters const& params) const { - if constexpr (std::is_same>::value) { // MT must be in RowMajor format - 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(); - } - } 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, void> { // MT must be in RowMajor format + 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(); } } diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index c593b9e6f..604ff202e 100644 --- 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 DenseMT; 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 @@ -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++; } } From 146c5d7efa6b722b6841a136d760d8e0cfdf90ac Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 7 Aug 2024 13:47:16 +0200 Subject: [PATCH 06/15] minor optimization --- include/convex_bodies/hpolytope.h | 5 +++-- include/random_walks/uniform_accelerated_billiard_walk.hpp | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index f27e7c934..99036ef0c 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -570,9 +570,9 @@ class HPolytope { template std::pair line_positive_intersect(Point const& r, - Point const& v, VT& Ar, VT& Av, + NT const& lambda_prev, VT& distances_vec, set_type& distances_set, AA_type const& AA, @@ -593,7 +593,8 @@ class HPolytope { *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); *(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value(); - distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); + if(*(dvec_data + it.row()) > params.moved_dist - lambda_prev) + distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); *(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row()); diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 604ff202e..76cf4c27f 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -159,7 +159,7 @@ struct AcceleratedBilliardWalk { std::pair pbpair; if constexpr (std::is_same>::value) { - pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _distances_vec, + pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, _distances_vec, _distances_set, _AA, _update_parameters); } else { pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev, From 76f1a63484147af7e3d19c42bae48257179e6bac Mon Sep 17 00:00:00 2001 From: lucaperju Date: Thu, 8 Aug 2024 13:51:52 +0200 Subject: [PATCH 07/15] replace std set with heap --- include/convex_bodies/hpolytope.h | 41 +++++-- .../uniform_accelerated_billiard_walk.hpp | 114 +++++++++++++++++- 2 files changed, 138 insertions(+), 17 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 99036ef0c..bd13b4a97 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -585,31 +585,50 @@ class HPolytope { // real distances_vec = distances_vec - params.moved_dist NT* Av_data = Av.data(); - NT* Ar_data = Ar.data(); - NT* dvec_data = distances_vec.data(); + //NT* Ar_data = Ar.data(); + //NT* dvec_data = distances_vec.data(); const NT* b_data = b.data(); + for (Eigen::SparseMatrix::InnerIterator it(AA, params.facet_prev); it; ++it) { + // NT old_Av = *(Av_data + it.row()); + *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); - *(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value(); + //*(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value(); - if(*(dvec_data + it.row()) > params.moved_dist - lambda_prev) - distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); + //if(*(dvec_data + it.row()) > params.moved_dist - lambda_prev) + // distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); - *(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row()); + // *(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row()); + // *(dvec_data + it.row()) = (*(dvec_data + it.row()) - params.moved_dist) * old_Av / *(Av_data + it.row()) + params.moved_dist; + // *(dvec_data + it.row()) += (*(dvec_data + it.row()) - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); + + NT val = distances_set.get_val(it.row()); + val += (val - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); + distances_set.change_val(it.row(), val, params.moved_dist); - if(*(dvec_data + it.row()) > params.moved_dist) - distances_set.insert(std::make_pair(*(dvec_data + it.row()), it.row())); + //if(*(dvec_data + it.row()) > params.moved_dist) + // distances_set.insert(std::make_pair(*(dvec_data + it.row()), it.row())); } + /* + std::cout << "got here post!!" << std::endl; + std::cout << std::endl; + std::cout << distances_set.heap_size << std::endl; + std::cout << params.moved_dist << std::endl; + std::cout << lambda_prev << std::endl; + for(int i = 0; i < num_of_hyperplanes(); ++i) { + std::cout << distances_set.heap[i].first << ' '; + } + std::cout << std::endl << std::endl;*/ - auto it = distances_set.upper_bound(std::make_pair(params.moved_dist, 0)); + /*auto it = distances_set.upper_bound(std::make_pair(-9999999, 0)); if(it == distances_set.end()) { std::cout << "something went wrong when trying to get lowest positive value" << std::endl; throw "all values from the set were negative"; } - - std::pair ans = (*it); + */ + std::pair ans = distances_set.get_min(); ans.first -= params.moved_dist; params.inner_vi_ak = *(Av_data + ans.second); diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 76cf4c27f..075f771e5 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -14,6 +14,103 @@ #include "sampling/sphere.hpp" #include #include +#include + +template +class Heap { +public: + int n, heap_size; + std::vector> heap; + std::vector> vec; + +private: + void 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) { + child += 1; + } + if(heap[child].first < heap[index].first) + { + std::swap(heap[child], heap[index]); + std::swap(vec[heap[child].second].second, vec[heap[index].second].second); + index = child; + } else { + return; + } + } + } + + void siftUp(int index) { + while(index > 0 && heap[(index - 1) >> 1].first > 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; + } + } + +public: + Heap() {} + + Heap(int n) : n(n), heap_size(0) { + heap.resize(n); + vec.resize(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 > 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); + } + } + + NT get_val (const int &index) { + return vec[index].first; + } + + std::pair get_min () { + return heap[0]; + } + + void remove (const int index) { // takes the index from the heap + 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; + siftDown(index); + } + + void insert (const std::pair val) { + heap[heap_size++] = val; + siftUp(heap_size - 1); + } + + void change_val(const int& index, const NT& new_val, const NT& moved_dist) { // takes the index from the vec + if(new_val < moved_dist) { // should not be inserted into the heap + remove(vec[index].second); + } else { // should be inserted into the heap + if(vec[index].second == -1) { + vec[index].second = heap_size; + insert({new_val, index}); + } else { + heap[vec[index].second].first = new_val; + siftDown(vec[index].second); + siftUp(vec[index].second); + } + } + vec[index].first = new_val; + } +}; // Billiard walk which accelarates each step for uniform distribution @@ -136,18 +233,21 @@ struct AcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; if constexpr (std::is_same>::value) { _update_parameters.moved_dist = _lambda_prev; - _distances_set.clear(); + //_distances_set.clear(); _distances_vec.setZero(P.num_of_hyperplanes()); typename Point::Coeff b = P.get_vec(); NT* b_data = b.data(); - NT* dvec_data = _distances_vec.data(); + //NT* dvec_data = _distances_vec.data(); NT* Ar_data = _lambdas.data(); NT* Av_data = _Av.data(); for(int i = 0; i < P.num_of_hyperplanes(); ++i) { - *(dvec_data + i) = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); - if(*(dvec_data + i) > _update_parameters.moved_dist) - _distances_set.insert(std::make_pair(*(dvec_data + i), i)); + _distances_set.vec[i].first = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + //std::cout << _distances_set.vec[i].first << ' '; + //*(dvec_data + i) = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); + //if(*(dvec_data + i) > _update_parameters.moved_dist) + // _distances_set.insert(std::make_pair(*(dvec_data + i), i)); } + _distances_set.rebuild(_update_parameters.moved_dist); } else { _p += (_lambda_prev * _v); } @@ -283,6 +383,7 @@ struct AcceleratedBilliardWalk _Av.setZero(P.num_of_hyperplanes()); _p = p; _v = GetDirection::apply(n, rng); + _distances_set = Heap(P.num_of_hyperplanes()); NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; @@ -350,7 +451,8 @@ struct AcceleratedBilliardWalk typename Point::Coeff _lambdas; typename Point::Coeff _Av; typename Point::Coeff _distances_vec; - std::set> _distances_set; + Heap _distances_set; + //std::set> _distances_set; }; }; From 67b0e0edd8267d840a94befcf19a3595b31c0b9c Mon Sep 17 00:00:00 2001 From: lucaperju Date: Thu, 8 Aug 2024 13:55:08 +0200 Subject: [PATCH 08/15] remove comments --- include/convex_bodies/hpolytope.h | 40 +------------------ .../uniform_accelerated_billiard_walk.hpp | 11 +---- 2 files changed, 2 insertions(+), 49 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index bd13b4a97..979703137 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -573,61 +573,23 @@ class HPolytope { VT& Ar, VT& Av, NT const& lambda_prev, - VT& distances_vec, set_type& distances_set, AA_type const& AA, update_parameters& params) const { NT inner_prev = params.inner_vi_ak; - - // real Ar = Ar + params.moved_dist * Av - // real r = r + params.moved_dist * v - // real distances_vec = distances_vec - params.moved_dist NT* Av_data = Av.data(); - //NT* Ar_data = Ar.data(); - //NT* dvec_data = distances_vec.data(); const NT* b_data = b.data(); for (Eigen::SparseMatrix::InnerIterator it(AA, params.facet_prev); it; ++it) { - // NT old_Av = *(Av_data + it.row()); - *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); - //*(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value(); - - //if(*(dvec_data + it.row()) > params.moved_dist - lambda_prev) - // distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row())); - - // *(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row()); - // *(dvec_data + it.row()) = (*(dvec_data + it.row()) - params.moved_dist) * old_Av / *(Av_data + it.row()) + params.moved_dist; - // *(dvec_data + it.row()) += (*(dvec_data + it.row()) - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); - NT val = distances_set.get_val(it.row()); val += (val - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); distances_set.change_val(it.row(), val, params.moved_dist); - - //if(*(dvec_data + it.row()) > params.moved_dist) - // distances_set.insert(std::make_pair(*(dvec_data + it.row()), it.row())); } - /* - std::cout << "got here post!!" << std::endl; - std::cout << std::endl; - std::cout << distances_set.heap_size << std::endl; - std::cout << params.moved_dist << std::endl; - std::cout << lambda_prev << std::endl; - for(int i = 0; i < num_of_hyperplanes(); ++i) { - std::cout << distances_set.heap[i].first << ' '; - } - std::cout << std::endl << std::endl;*/ - - /*auto it = distances_set.upper_bound(std::make_pair(-9999999, 0)); - - if(it == distances_set.end()) { - std::cout << "something went wrong when trying to get lowest positive value" << std::endl; - throw "all values from the set were negative"; - } - */ + std::pair ans = distances_set.get_min(); ans.first -= params.moved_dist; diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 075f771e5..960871017 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -233,19 +233,12 @@ struct AcceleratedBilliardWalk _lambda_prev = dl * pbpair.first; if constexpr (std::is_same>::value) { _update_parameters.moved_dist = _lambda_prev; - //_distances_set.clear(); - _distances_vec.setZero(P.num_of_hyperplanes()); typename Point::Coeff b = P.get_vec(); NT* b_data = b.data(); - //NT* dvec_data = _distances_vec.data(); 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)); - //std::cout << _distances_set.vec[i].first << ' '; - //*(dvec_data + i) = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i)); - //if(*(dvec_data + i) > _update_parameters.moved_dist) - // _distances_set.insert(std::make_pair(*(dvec_data + i), i)); } _distances_set.rebuild(_update_parameters.moved_dist); } else { @@ -259,7 +252,7 @@ struct AcceleratedBilliardWalk { std::pair pbpair; if constexpr (std::is_same>::value) { - pbpair = P.line_positive_intersect(_p, _lambdas, _Av, _lambda_prev, _distances_vec, + 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, @@ -450,9 +443,7 @@ struct AcceleratedBilliardWalk update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - typename Point::Coeff _distances_vec; Heap _distances_set; - //std::set> _distances_set; }; }; From 0f76307c350d4d42fa574575f412fabcc38affb7 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 9 Aug 2024 11:42:50 +0200 Subject: [PATCH 09/15] fix bug --- .../uniform_accelerated_billiard_walk.hpp | 21 ++++++++++++------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 960871017..e3a63fe06 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -16,6 +16,8 @@ #include #include +const double eps = 0.000000001; + template class Heap { public: @@ -24,29 +26,31 @@ class Heap { std::vector> vec; private: - void siftDown(int index) { + 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) { + if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) { child += 1; } - if(heap[child].first < heap[index].first) + 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; + return index; } } + return index; } - void siftUp(int index) { - while(index > 0 && heap[(index - 1) >> 1].first > heap[index].first) { + 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; } public: @@ -79,7 +83,7 @@ class Heap { return heap[0]; } - void remove (const int index) { // takes the index from the heap + void remove (int index) { // takes the index from the heap if(index == -1) { return; } @@ -87,7 +91,8 @@ class Heap { 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; - siftDown(index); + index = siftDown(index); + siftUp(index); } void insert (const std::pair val) { From fe5641df49dd5fe205f26e3c39ffbc98066307ae Mon Sep 17 00:00:00 2001 From: lucaperju Date: Wed, 14 Aug 2024 12:29:17 +0200 Subject: [PATCH 10/15] minor changes --- .../uniform_accelerated_billiard_walk.hpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index e3a63fe06..d1d7fc056 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -65,7 +65,7 @@ class Heap { heap_size = 0; for(int i = 0; i < n; ++i) { vec[i].second = -1; - if(vec[i].first > moved_dist) { + if(vec[i].first - eps > moved_dist) { vec[i].second = heap_size; heap[heap_size++] = {vec[i].first, i}; } @@ -83,7 +83,8 @@ class Heap { return heap[0]; } - void remove (int index) { // takes the index from the heap + void remove (int index) { // takes the index from the vec + index = vec[index].second; if(index == -1) { return; } @@ -96,24 +97,26 @@ class Heap { } 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); } void change_val(const int& index, const NT& new_val, const NT& moved_dist) { // takes the index from the vec - if(new_val < moved_dist) { // should not be inserted into the heap - remove(vec[index].second); - } else { // should be inserted into the heap + if(new_val < moved_dist - eps) { + vec[index].first = new_val; + remove(index); + } else { if(vec[index].second == -1) { - vec[index].second = heap_size; 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); } } - vec[index].first = new_val; } }; From d9bd3c2268fdcce2854b8715ce3c59357ccab9b9 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 16 Aug 2024 09:48:20 +0200 Subject: [PATCH 11/15] minor changes and added a few explanatory comments --- include/convex_bodies/hpolytope.h | 32 ++++++--- .../uniform_accelerated_billiard_walk.hpp | 66 +++++++++++-------- 2 files changed, 63 insertions(+), 35 deletions(-) diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 979703137..9cabcd9ab 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -578,29 +578,41 @@ class HPolytope { update_parameters& params) const { NT inner_prev = params.inner_vi_ak; - NT* Av_data = Av.data(); - const NT* b_data = b.data(); + + // Av += (-2.0 * inner_prev) * AA.col(params.facet_prev) for (Eigen::SparseMatrix::InnerIterator it(AA, params.facet_prev); it; ++it) { + + // 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) + *(Av_data + it.row()) += (-2.0 * inner_prev) * it.value(); + + // 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 + + // new_val(row) = (old_val(row) - params.moved_dist) * old_Av(row) / new_Av(row) + params.moved_dist; + // new_val(row) = (old_val(row) - params.moved_dist) * (new_Av(row) + 2.0 * inner_prev * it.value()) / new_Av(row) + params.moved_dist; + // new_val(row) = (old_val(row) - params.moved_dist) * (1 + (2.0 * inner_prev * it.value()) / new_Av(row) ) + params.moved_dist; + + // new_val(row) = old_val(row) + (old_val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / new_Av(row) + + // val(row) += (val(row) - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); + NT val = distances_set.get_val(it.row()); val += (val - params.moved_dist) * 2.0 * inner_prev * it.value() / *(Av_data + it.row()); distances_set.change_val(it.row(), val, params.moved_dist); } - + std::pair ans = distances_set.get_min(); ans.first -= params.moved_dist; - params.inner_vi_ak = *(Av_data + ans.second); params.facet_prev = ans.second; - /*if(ans.first < 0.00000001) { - std::cout << "distance of 0 found" << std::endl; - exit(0); - }*/ - return ans; } @@ -996,6 +1008,8 @@ class HPolytope { 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 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 diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index d1d7fc056..97429b308 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -16,10 +16,13 @@ #include #include -const double eps = 0.000000001; +const double eps = 0.0000000001; +// 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 Heap { +class BoundaryOracleHeap { public: int n, heap_size; std::vector> heap; @@ -53,6 +56,28 @@ class Heap { 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: Heap() {} @@ -61,6 +86,8 @@ class Heap { 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) { @@ -75,35 +102,21 @@ class Heap { } } + // 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]; } - void remove (int index) { // takes the index from the vec - 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); - } - - 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); - } - - void change_val(const int& index, const NT& new_val, const NT& moved_dist) { // takes the index from the vec + // 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); @@ -222,6 +235,8 @@ struct AcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; + typename Point::Coeff b = P.get_vec(); + NT* b_data = b.data(); for (auto j=0u; j>::value) { _update_parameters.moved_dist = _lambda_prev; - typename Point::Coeff b = P.get_vec(); - NT* b_data = b.data(); 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); @@ -451,7 +465,7 @@ struct AcceleratedBilliardWalk update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - Heap _distances_set; + BoundaryOracleHeap _distances_set; }; }; From 3189a0b91cd85e8d315e0bd809a8f7b424be6185 Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 16 Aug 2024 10:23:02 +0200 Subject: [PATCH 12/15] fix bug --- include/random_walks/uniform_accelerated_billiard_walk.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 97429b308..b27dd4b81 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -79,9 +79,9 @@ class BoundaryOracleHeap { } public: - Heap() {} + BoundaryOracleHeap() {} - Heap(int n) : n(n), heap_size(0) { + BoundaryOracleHeap(int n) : n(n), heap_size(0) { heap.resize(n); vec.resize(n); } @@ -398,7 +398,7 @@ struct AcceleratedBilliardWalk _Av.setZero(P.num_of_hyperplanes()); _p = p; _v = GetDirection::apply(n, rng); - _distances_set = Heap(P.num_of_hyperplanes()); + _distances_set = BoundaryOracleHeap(P.num_of_hyperplanes()); NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; From 5e4990434d131e4220ea0e9d36d9ffb4ceeeee3d Mon Sep 17 00:00:00 2001 From: lucaperju Date: Fri, 16 Aug 2024 10:44:13 +0200 Subject: [PATCH 13/15] fix bug --- include/random_walks/uniform_accelerated_billiard_walk.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index b27dd4b81..30f85e584 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -235,8 +235,10 @@ struct AcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; - typename Point::Coeff b = P.get_vec(); - NT* b_data = b.data(); + if constexpr (std::is_same>::value) { + typename Point::Coeff b = P.get_vec(); + NT* b_data = b.data(); + } for (auto j=0u; j Date: Fri, 16 Aug 2024 10:50:10 +0200 Subject: [PATCH 14/15] fix b_data bug --- include/random_walks/uniform_accelerated_billiard_walk.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 30f85e584..d9bec3f6a 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -235,9 +235,11 @@ struct AcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; + typename Point::Coeff b; + NT* b_data; if constexpr (std::is_same>::value) { - typename Point::Coeff b = P.get_vec(); - NT* b_data = b.data(); + b = P.get_vec(); + b_data = b.data(); } for (auto j=0u; j Date: Mon, 19 Aug 2024 11:11:53 +0200 Subject: [PATCH 15/15] change epsilon notation --- include/random_walks/uniform_accelerated_billiard_walk.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index d9bec3f6a..f91abee4b 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -16,7 +16,7 @@ #include #include -const double eps = 0.0000000001; +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