Skip to content

Commit

Permalink
minor changes and added a few explanatory comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucaperju
lucaperju committed Aug 16, 2024
1 parent fe5641d commit d9bd3c2
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 35 deletions.
32 changes: 23 additions & 9 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>::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<NT, int> 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;
}

Expand Down Expand Up @@ -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 <typename update_parameters>
auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
-> std::enable_if_t<std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>> && !std::is_same_v<update_parameters, int>, void> { // MT must be in RowMajor format
Expand Down
66 changes: 40 additions & 26 deletions include/random_walks/uniform_accelerated_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,13 @@
#include <set>
#include <vector>

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<typename NT>
class Heap {
class BoundaryOracleHeap {
public:
int n, heap_size;
std::vector<std::pair<NT, int>> heap;
Expand Down Expand Up @@ -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<NT, int> val) {
vec[val.second].second = heap_size;
vec[val.second].first = val.first;
heap[heap_size++] = val;
siftUp(heap_size - 1);
}

public:
Heap() {}

Expand All @@ -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) {
Expand All @@ -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<NT, int> 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<NT, int> 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);
Expand Down Expand Up @@ -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<walk_length; ++j)
{
Expand All @@ -241,13 +256,12 @@ struct AcceleratedBilliardWalk
_lambda_prev = dl * pbpair.first;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::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);
Expand Down Expand Up @@ -451,7 +465,7 @@ struct AcceleratedBilliardWalk
update_parameters _update_parameters;
typename Point::Coeff _lambdas;
typename Point::Coeff _Av;
Heap<NT> _distances_set;
BoundaryOracleHeap<NT> _distances_set;
};

};
Expand Down

0 comments on commit d9bd3c2

Please sign in to comment.