Skip to content

Commit

Permalink
finish sparse optimizations for UniformABW
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucaperju
lucaperju committed Aug 7, 2024
1 parent 37f0ca0 commit 725606b
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 22 deletions.
60 changes: 58 additions & 2 deletions include/convex_bodies/hpolytope.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -525,7 +525,7 @@ class HPolytope {
VT& Av,
NT const& lambda_prev,
DenseMT const& AA,
update_parameters& params) const//
update_parameters& params) const
{

NT min_plus = std::numeric_limits<NT>::max();
Expand Down Expand Up @@ -567,6 +567,62 @@ class HPolytope {
}



template <typename update_parameters, typename set_type, typename AA_type>
std::pair<NT, int> line_positive_intersect(Point const& r,
Point const& v,
VT& Ar,
VT& Av,
VT& distances_vec,
set_type& distances_set,
AA_type const& AA,
update_parameters& params) const
{
NT inner_prev = params.inner_vi_ak;

// real Ar = Ar + params.moved_dist * Av
// real r = r + params.moved_dist * v
// real distances_vec = distances_vec - params.moved_dist

NT* Av_data = Av.data();
NT* Ar_data = Ar.data();
NT* dvec_data = distances_vec.data();
const NT* b_data = b.data();
for (Eigen::SparseMatrix<double>::InnerIterator it(AA, params.facet_prev); it; ++it) {

*(Av_data + it.row()) += (-2.0 * inner_prev) * it.value();
*(Ar_data + it.row()) -= (-2.0 * inner_prev * params.moved_dist) * it.value();

distances_set.erase(std::make_pair(*(dvec_data + it.row()), it.row()));

*(dvec_data + it.row()) = (*(b_data + it.row()) - *(Ar_data + it.row())) / *(Av_data + it.row());

if(*(dvec_data + it.row()) > params.moved_dist)
distances_set.insert(std::make_pair(*(dvec_data + it.row()), it.row()));
}

auto it = distances_set.upper_bound(std::make_pair(params.moved_dist, 0));

if(it == distances_set.end()) {
std::cout << "something went wrong when trying to get lowest positive value" << std::endl;
throw "all values from the set were negative";
}

std::pair<NT, int> ans = (*it);
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;
}


template <typename update_parameters>
std::pair<NT, int> line_positive_intersect(Point const& r,
Point const& v,
Expand Down Expand Up @@ -954,7 +1010,7 @@ class HPolytope {

template <typename update_parameters>
void compute_reflection(Point &v, Point &p, update_parameters const& params) const {
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) { // is faster only if MT is in RowMajor format
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) { // MT must be in RowMajor format
NT* v_data = v.pointerToData();
NT* p_data = p.pointerToData();
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
Expand Down
63 changes: 43 additions & 20 deletions include/random_walks/uniform_accelerated_billiard_walk.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@

#include "sampling/sphere.hpp"
#include <Eigen/Eigen>
#include <set>


// Billiard walk which accelarates each step for uniform distribution
Expand Down Expand Up @@ -62,6 +63,8 @@ struct AcceleratedBilliardWalk
typedef typename Polytope::MT MT;
typedef typename Polytope::DenseMT DenseMT;
typedef typename Point::FT NT;
using AA_type = std::conditional_t< std::is_same_v<MT, typename Eigen::SparseMatrix<NT, Eigen::RowMajor>>, typename Eigen::SparseMatrix<NT>, DenseMT >;
// AA is sparse colMajor if MT is sparse rowMajor, and Dense otherwise

template <typename GenericPolytope>
Walk(GenericPolytope &P, Point const& p, RandomNumberGenerator &rng)
Expand All @@ -72,7 +75,11 @@ struct AcceleratedBilliardWalk
_update_parameters = update_parameters();
_L = compute_diameter<GenericPolytope>
::template compute<NT>(P);
if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::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);
}
Expand All @@ -88,7 +95,11 @@ struct AcceleratedBilliardWalk
_L = params.set_L ? params.m_L
: compute_diameter<GenericPolytope>
::template compute<NT>(P);
if constexpr (std::is_same<AA_type, Eigen::SparseMatrix<NT>>::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);
}
Expand All @@ -114,13 +125,7 @@ struct AcceleratedBilliardWalk
Point p0 = _p;

it = 0;
std::pair<NT, int> 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<NT, int> pbpair = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);

if (T <= pbpair.first) {
_p += (T * _v);
Expand All @@ -129,18 +134,37 @@ struct AcceleratedBilliardWalk
}

_lambda_prev = dl * pbpair.first;
_p += (_lambda_prev * _v);
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
_update_parameters.moved_dist = _lambda_prev;
_distances_set.clear();
_distances_vec.setZero(P.num_of_hyperplanes());
typename Point::Coeff b = P.get_vec();
NT* b_data = b.data();
NT* dvec_data = _distances_vec.data();
NT* Ar_data = _lambdas.data();
NT* Av_data = _Av.data();
for(int i = 0; i < P.num_of_hyperplanes(); ++i) {
*(dvec_data + i) = ( *(b_data + i) - (*(Ar_data + i)) ) / (*(Av_data + i));
if(*(dvec_data + i) > _update_parameters.moved_dist)
_distances_set.insert(std::make_pair(*(dvec_data + i), i));
}
} else {
_p += (_lambda_prev * _v);
}
T -= _lambda_prev;
P.compute_reflection(_v, _p, _update_parameters);
it++;

_update_parameters.moved_dist = 0.0;

while (it < _rho)
{
std::pair<NT, int> pbpair
= P.line_positive_intersect(_p, _v, _lambdas, _Av, _lambda_prev,
_AA, _update_parameters); //
{
std::pair<NT, int> pbpair;
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
pbpair = P.line_positive_intersect(_p, _v, _lambdas, _Av, _distances_vec,
_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;
Expand All @@ -150,17 +174,16 @@ struct AcceleratedBilliardWalk
if constexpr (std::is_same<MT, Eigen::SparseMatrix<NT, Eigen::RowMajor>>::value) {
_update_parameters.moved_dist += _lambda_prev;
} else {
_p += (_lambda_prev * _v);// done
_p += (_lambda_prev * _v);
}
T -= _lambda_prev;
P.compute_reflection(_v, _p, _update_parameters);// done
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;
Expand Down Expand Up @@ -256,7 +279,6 @@ 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;
Expand Down Expand Up @@ -322,12 +344,13 @@ struct AcceleratedBilliardWalk
Point _p;
Point _v;
NT _lambda_prev;
DenseMT _AA;
AA_type _AA;
unsigned int _rho;
update_parameters _update_parameters;
typename Point::Coeff _lambdas;
typename Point::Coeff _Av;
bool was_reset;
typename Point::Coeff _distances_vec;
std::set<std::pair<NT, int>> _distances_set;
};

};
Expand Down

0 comments on commit 725606b

Please sign in to comment.