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; }; };