Skip to content

Commit 4fc2b3f

Browse files
lucaperjulucaperju
committed
minor changes
1 parent 19565dc commit 4fc2b3f

File tree

4 files changed

+148
-135
lines changed

4 files changed

+148
-135
lines changed

‎include/convex_bodies/hpolytope.h

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1024,8 +1024,8 @@ class HPolytope {
10241024

10251025
// function to compute reflection for GaBW random walk
10261026
// compatible when the polytope is both dense or sparse
1027-
template <typename AE_type, typename update_parameters>
1028-
NT compute_reflection(Point &v, Point &p, AE_type const &AE, VT const &AEA, NT &vEv, update_parameters const &params) const {
1027+
template <typename DenseSparseMT, typename update_parameters>
1028+
NT compute_reflection(Point &v, Point &p, NT &vEv, DenseSparseMT const &AE, VT const &AEA, update_parameters const &params) const {
10291029

10301030
NT new_vEv;
10311031
if constexpr (!std::is_same_v<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
@@ -1035,7 +1035,7 @@ class HPolytope {
10351035
v += a;
10361036
} else {
10371037

1038-
if constexpr(!std::is_same_v<AE_type, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
1038+
if constexpr(!std::is_same_v<DenseSparseMT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>) {
10391039
VT x = v.getCoefficients();
10401040
new_vEv = vEv - (4.0 * params.inner_vi_ak) * (AE.row(params.facet_prev).dot(x) - params.inner_vi_ak * AEA(params.facet_prev));
10411041
} else {

‎include/random_walks/accelerated_billiard_walk_utils.hpp

Lines changed: 133 additions & 0 deletions
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

‎include/random_walks/gaussian_accelerated_billiard_walk.hpp

Lines changed: 11 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
#include "generators/boost_random_number_generator.hpp"
2222
#include "sampling/ellipsoid.hpp"
2323
#include "random_walks/uniform_billiard_walk.hpp"
24+
#include "random_walks/accelerated_billiard_walk_utils.hpp"
2425

2526
#include "random_walks/compute_diameter.hpp"
2627

@@ -72,11 +73,10 @@ struct GaussianAcceleratedBilliardWalk
7273
typedef typename Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> DenseMT;
7374
typedef typename Polytope::VT VT;
7475
typedef typename Point::FT NT;
75-
using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >;
7676
// AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise
77-
using AE_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>> && std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type>::value, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>, DenseMT >;
78-
// ( ( ( )) ( ( ) ) ( ) )
77+
using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >;
7978
// AE is sparse rowMajor if (MT is sparse rowMajor and E is sparse), and Dense otherwise
79+
using AE_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>> && std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type>::value, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>, DenseMT >;
8080

8181
void computeLcov(const E_type E)
8282
{
@@ -163,12 +163,6 @@ struct GaussianAcceleratedBilliardWalk
163163
int it;
164164
NT coef = 1.0;
165165
NT vEv;
166-
typename Point::Coeff b;
167-
NT* b_data;
168-
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
169-
b = P.get_vec();
170-
b_data = b.data();
171-
}
172166

173167
for (auto j=0u; j<walk_length; ++j)
174168
{
@@ -193,6 +187,10 @@ struct GaussianAcceleratedBilliardWalk
193187

194188
_lambda_prev = dl * pbpair.first;
195189
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
190+
typename Point::Coeff b;
191+
NT* b_data;
192+
b = P.get_vec();
193+
b_data = b.data();
196194
_update_parameters.moved_dist = _lambda_prev;
197195
NT* Ar_data = _lambdas.data();
198196
NT* Av_data = _Av.data();
@@ -207,7 +205,7 @@ struct GaussianAcceleratedBilliardWalk
207205
T -= _lambda_prev;
208206

209207
T = T * coef;
210-
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
208+
coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
211209
T = T / coef;
212210

213211
it++;
@@ -237,7 +235,7 @@ struct GaussianAcceleratedBilliardWalk
237235
T -= _lambda_prev;
238236

239237
T = T * coef;
240-
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
238+
coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
241239
T = T / coef;
242240

243241
it++;
@@ -298,7 +296,7 @@ struct GaussianAcceleratedBilliardWalk
298296
T -= _lambda_prev;
299297

300298
T = T * coef;
301-
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
299+
coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
302300
T = T / coef;
303301

304302
while (it <= _rho)
@@ -320,7 +318,7 @@ struct GaussianAcceleratedBilliardWalk
320318
T -= _lambda_prev;
321319

322320
T = T * coef;
323-
coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
321+
coef = P.compute_reflection(_v, _p, vEv, _AE, _AEA, _update_parameters);
324322
T = T / coef;
325323

326324
it++;

‎include/random_walks/uniform_accelerated_billiard_walk.hpp

Lines changed: 1 addition & 119 deletions
Original file line numberDiff line numberDiff line change
@@ -13,125 +13,7 @@
1313

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

13618

13719
// Billiard walk which accelarates each step for uniform distribution

0 commit comments

Comments
 (0)