Skip to content

Commit 5007b1d

Browse files
lucaperjulucaperju
authored
Sparse optimizations for GaBW sampling (#331)
* simplify gabw sampling * optimize gabw for sparse polytopes * minor changes * add SPARSE boolean for clarity * compute_reflection_abw_sparse and improve comments
1 parent 9cb7226 commit 5007b1d

6 files changed

+322
-190
lines changed

‎.gitignore

+2-1
Original file line numberDiff line numberDiff line change
@@ -9,4 +9,5 @@ test/Testing/Temporary/CTestCostData.txt
99
build/
1010
.vscode
1111
.DS_Store
12-
.cache
12+
.cache
13+
shell.nix

‎include/convex_bodies/hpolytope.h

+52-17
Original file line numberDiff line numberDiff line change
@@ -580,17 +580,23 @@ class HPolytope {
580580
NT inner_prev = params.inner_vi_ak;
581581
NT* Av_data = Av.data();
582582

583+
// Updating Av due to the change in direction caused by the previous reflection
583584
// Av += (-2.0 * inner_prev) * AA.col(params.facet_prev)
584585

585586
for (Eigen::SparseMatrix<double>::InnerIterator it(AA, params.facet_prev); it; ++it) {
586587

588+
// val(row) - params.moved_dist = The distance until we would hit the facet given by row
589+
// all those values are stored inside distances_set for quick retrieval of the minimum
587590

591+
// Before updating Av(row)
588592
// val(row) = (b(row) - Ar(row)) / Av(row) + params.moved_dist
589593
// (val(row) - params.moved_dist) = (b(row) - Ar(row)) / Av(row)
590594
// (val(row) - params.moved_dist) * Av(row) = b(row) - Ar(row)
595+
// b(row) - Ar(row) = (val(row) - params.moved_dist) * Av(row)
591596

592597
*(Av_data + it.row()) += (-2.0 * inner_prev) * it.value();
593598

599+
// After updating Av(row)
594600
// b(row) - Ar(row) = (old_val(row) - params.moved_dist) * old_Av(row)
595601
// new_val(row) = (b(row) - Ar(row) ) / new_Av(row) + params.moved_dist
596602
// new_val(row) = ((old_val(row) - params.moved_dist) * old_Av(row)) / new_Av(row) + params.moved_dist
@@ -608,8 +614,12 @@ class HPolytope {
608614
distances_set.change_val(it.row(), val, params.moved_dist);
609615
}
610616

617+
// Finding the distance to the closest facet and its row
611618
std::pair<NT, int> ans = distances_set.get_min();
619+
620+
// Subtract params.moved_dist to obtain the actual distance to the facet
612621
ans.first -= params.moved_dist;
622+
613623
params.inner_vi_ak = *(Av_data + ans.second);
614624
params.facet_prev = ans.second;
615625

@@ -1002,34 +1012,59 @@ class HPolytope {
10021012
return total;
10031013
}
10041014

1015+
// Updates the velocity vector v and the position vector p after a reflection
10051016
template <typename update_parameters>
10061017
void compute_reflection(Point &v, Point const&, update_parameters const& params) const {
10071018
Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
10081019
v += a;
10091020
}
10101021

1011-
// updates the velocity vector v and the position vector p after a reflection
1012-
// the real value of p is given by p + moved_dist * v
1022+
// Only to be called when MT is in RowMajor format
1023+
// The real value of p is given by p + params.moved_dist * v
10131024
template <typename update_parameters>
1014-
auto compute_reflection(Point &v, Point &p, update_parameters const& params) const
1015-
-> 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
1016-
NT* v_data = v.pointerToData();
1017-
NT* p_data = p.pointerToData();
1018-
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
1019-
*(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
1020-
*(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
1021-
}
1025+
auto compute_reflection_abw_sparse(Point &v, Point &p, update_parameters const& params) const
1026+
{
1027+
NT* v_data = v.pointerToData();
1028+
NT* p_data = p.pointerToData();
1029+
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
1030+
*(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
1031+
*(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
1032+
}
10221033
}
10231034

1024-
template <typename update_parameters>
1025-
NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) const {
1035+
// function to compute reflection for GaBW random walk
1036+
// compatible when the polytope is both dense or sparse
1037+
template <typename DenseSparseMT, typename update_parameters>
1038+
NT compute_reflection(Point &v, Point &p, NT &vEv, DenseSparseMT const &AE, VT const &AEA, update_parameters const &params) const {
1039+
1040+
NT new_vEv;
1041+
if constexpr (!std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
1042+
Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
1043+
VT x = v.getCoefficients();
1044+
new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
1045+
v += a;
1046+
} else {
1047+
1048+
if constexpr(!std::is_same_v<DenseSparseMT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
1049+
VT x = v.getCoefficients();
1050+
new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
1051+
} else {
1052+
new_vEv = vEv + 4.0 * params.inner_vi_ak * params.inner_vi_ak * AEA(params.facet_prev);
1053+
NT* v_data = v.pointerToData();
1054+
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(AE, params.facet_prev); it; ++it) {
1055+
new_vEv -= 4.0 * params.inner_vi_ak * it.value() * *(v_data + it.col());
1056+
}
1057+
}
10261058

1027-
Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
1028-
VT x = v.getCoefficients();
1029-
NT new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
1030-
v += a;
1059+
NT* v_data = v.pointerToData();
1060+
NT* p_data = p.pointerToData();
1061+
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
1062+
*(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
1063+
*(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
1064+
}
1065+
}
10311066
NT coeff = std::sqrt(vEv / new_vEv);
1032-
v = v * coeff;
1067+
vEv = new_vEv;
10331068
return coeff;
10341069
}
10351070

‎include/random_walks/accelerated_billiard_walk_utils.hpp

+133
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
// VolEsti (volume computation and sampling library)
2+
3+
// Copyright (c) 2012-2020 Vissarion Fisikopoulos
4+
// Copyright (c) 2018-2020 Apostolos Chalkis
5+
// Copyright (c) 2024 Luca Perju
6+
7+
// Licensed under GNU LGPL.3, see LICENCE file
8+
9+
#ifndef ACCELERATED_BILLIARD_WALK_UTILS_HPP
10+
#define ACCELERATED_BILLIARD_WALK_UTILS_HPP
11+
12+
#include <Eigen/Eigen>
13+
#include <vector>
14+
15+
const double eps = 1e-10;
16+
17+
// data structure which maintains the values of (b - Ar)/Av, and can extract the minimum positive value and the facet associated with it
18+
// 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
19+
// The heap will only contain the values from vec which are greater than moved_dist (so that they are positive)
20+
template<typename NT>
21+
class BoundaryOracleHeap {
22+
public:
23+
int n, heap_size;
24+
std::vector<std::pair<NT, int>> heap;
25+
std::vector<std::pair<NT, int>> vec;
26+
27+
private:
28+
int siftDown(int index) {
29+
while((index << 1) + 1 < heap_size) {
30+
int child = (index << 1) + 1;
31+
if(child + 1 < heap_size && heap[child + 1].first < heap[child].first - eps) {
32+
child += 1;
33+
}
34+
if(heap[child].first < heap[index].first - eps)
35+
{
36+
std::swap(heap[child], heap[index]);
37+
std::swap(vec[heap[child].second].second, vec[heap[index].second].second);
38+
index = child;
39+
} else {
40+
return index;
41+
}
42+
}
43+
return index;
44+
}
45+
46+
int siftUp(int index) {
47+
while(index > 0 && heap[(index - 1) >> 1].first - eps > heap[index].first) {
48+
std::swap(heap[(index - 1) >> 1], heap[index]);
49+
std::swap(vec[heap[(index - 1) >> 1].second].second, vec[heap[index].second].second);
50+
index = (index - 1) >> 1;
51+
}
52+
return index;
53+
}
54+
55+
// takes the index of a facet, and (in case it is in the heap) removes it from the heap.
56+
void remove (int index) {
57+
index = vec[index].second;
58+
if(index == -1) {
59+
return;
60+
}
61+
std::swap(heap[heap_size - 1], heap[index]);
62+
std::swap(vec[heap[heap_size - 1].second].second, vec[heap[index].second].second);
63+
vec[heap[heap_size - 1].second].second = -1;
64+
heap_size -= 1;
65+
index = siftDown(index);
66+
siftUp(index);
67+
}
68+
69+
// inserts a new value into the heap, with its associated facet
70+
void insert (const std::pair<NT, int> val) {
71+
vec[val.second].second = heap_size;
72+
vec[val.second].first = val.first;
73+
heap[heap_size++] = val;
74+
siftUp(heap_size - 1);
75+
}
76+
77+
public:
78+
BoundaryOracleHeap() {}
79+
80+
BoundaryOracleHeap(int n) : n(n), heap_size(0) {
81+
heap.resize(n);
82+
vec.resize(n);
83+
}
84+
85+
// rebuilds the heap with the existing values from vec
86+
// O(n)
87+
void rebuild (const NT &moved_dist) {
88+
heap_size = 0;
89+
for(int i = 0; i < n; ++i) {
90+
vec[i].second = -1;
91+
if(vec[i].first - eps > moved_dist) {
92+
vec[i].second = heap_size;
93+
heap[heap_size++] = {vec[i].first, i};
94+
}
95+
}
96+
for(int i = heap_size - 1; i >= 0; --i) {
97+
siftDown(i);
98+
}
99+
}
100+
101+
// returns (b(i) - Ar(i))/Av(i) + moved_dist
102+
// O(1)
103+
NT get_val (const int &index) {
104+
return vec[index].first;
105+
}
106+
107+
// returns the nearest facet
108+
// O(1)
109+
std::pair<NT, int> get_min () {
110+
return heap[0];
111+
}
112+
113+
// changes the stored value for a given facet, and updates the heap accordingly
114+
// O(logn)
115+
void change_val(const int& index, const NT& new_val, const NT& moved_dist) {
116+
if(new_val < moved_dist - eps) {
117+
vec[index].first = new_val;
118+
remove(index);
119+
} else {
120+
if(vec[index].second == -1) {
121+
insert({new_val, index});
122+
} else {
123+
heap[vec[index].second].first = new_val;
124+
vec[index].first = new_val;
125+
siftDown(vec[index].second);
126+
siftUp(vec[index].second);
127+
}
128+
}
129+
}
130+
};
131+
132+
133+
#endif

0 commit comments

Comments
 (0)