diff --git a/include/convex_bodies/hpolytope.h b/include/convex_bodies/hpolytope.h index 5f37be325..9cabcd9ab 100644 --- a/include/convex_bodies/hpolytope.h +++ b/include/convex_bodies/hpolytope.h @@ -567,6 +567,56 @@ class HPolytope { } + + template + std::pair line_positive_intersect(Point const& r, + VT& Ar, + VT& Av, + NT const& lambda_prev, + set_type& distances_set, + AA_type const& AA, + update_parameters& params) const + { + NT inner_prev = params.inner_vi_ak; + NT* Av_data = Av.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; + + return ans; + } + + template std::pair line_positive_intersect(Point const& r, Point const& v, @@ -953,12 +1003,24 @@ class HPolytope { } template - void compute_reflection(Point &v, const Point &, 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; } + // 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 + 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(); + } + } + template NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const ¶ms) const { 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..f91abee4b 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -13,6 +13,125 @@ #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); + } + } + } +}; // Billiard walk which accelarates each step for uniform distribution @@ -39,12 +158,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,8 +178,11 @@ struct AcceleratedBilliardWalk struct Walk { typedef typename Polytope::PointType Point; - typedef typename Eigen::Matrix MT; + typedef typename Polytope::MT MT; + 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 template Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng) @@ -70,7 +193,11 @@ struct AcceleratedBilliardWalk _update_parameters = update_parameters(); _L = compute_diameter ::template compute(P); - _AA.noalias() = (MT)(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); } @@ -86,7 +213,11 @@ struct AcceleratedBilliardWalk _L = params.set_L ? params.m_L : compute_diameter ::template compute(P); - _AA.noalias() = (MT)(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); } @@ -104,6 +235,12 @@ struct AcceleratedBilliardWalk NT T; const NT dl = 0.995; int it; + 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); @@ -127,30 +258,51 @@ struct AcceleratedBilliardWalk } _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; P.compute_reflection(_v, _p, _update_parameters); it++; 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, _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); _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); + } T -= _lambda_prev; 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; @@ -246,11 +398,11 @@ 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; _v = GetDirection::apply(n, rng); + _distances_set = BoundaryOracleHeap(P.num_of_hyperplanes()); NT T = -std::log(rng.sample_urdist()) * _L; Point p0 = _p; @@ -312,12 +464,12 @@ struct AcceleratedBilliardWalk Point _p; Point _v; NT _lambda_prev; - MT _AA; + AA_type _AA; unsigned int _rho; update_parameters _update_parameters; typename Point::Coeff _lambdas; typename Point::Coeff _Av; - bool was_reset; + BoundaryOracleHeap _distances_set; }; }; 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;