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