Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sparse optimizations for UniformABW sampling #329

Merged
merged 15 commits into from
Aug 21, 2024
Merged

Sparse optimizations for UniformABW sampling #329

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
ae5b45f
support for eigen sparse rowmajor matrix
lucaperju Aug 5, 2024
37f0ca0
optimize updating the position in abw
lucaperju Aug 5, 2024
725606b
finish sparse optimizations for UniformABW
lucaperju Aug 7, 2024
c61c4bd
temporary fix to failing tests
lucaperju Aug 7, 2024
f066091
fix bug with reflection
lucaperju Aug 7, 2024
146c5d7
minor optimization
lucaperju Aug 7, 2024
76f1a63
replace std set with heap
lucaperju Aug 8, 2024
67b0e0e
remove comments
lucaperju Aug 8, 2024
0f76307
fix bug
lucaperju Aug 9, 2024
fe5641d
minor changes
lucaperju Aug 14, 2024
d9bd3c2
minor changes and added a few explanatory comments
lucaperju Aug 16, 2024
3189a0b
fix bug
lucaperju Aug 16, 2024
5e49904
fix bug
lucaperju Aug 16, 2024
b281839
fix b_data bug
lucaperju Aug 16, 2024
e141605
change epsilon notation
lucaperju Aug 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 64 additions & 2 deletions include/convex_bodies/hpolytope.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -567,6 +567,56 @@ class HPolytope {
}



template <typename update_parameters, typename set_type, typename AA_type>
std::pair<NT, int> line_positive_intersect(Point const& r,
VT& Ar,
VT& Av,
NT const& lambda_prev,
set_type& distances_set,
AA_type const& AA,
update_parameters& params) const
{
NT inner_prev = params.inner_vi_ak;
NT* Av_data = Av.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;

return ans;
}


template <typename update_parameters>
std::pair<NT, int> line_positive_intersect(Point const& r,
Point const& v,
Expand Down Expand Up @@ -953,12 +1003,24 @@ class HPolytope {
}

template <typename update_parameters>
void compute_reflection(Point &v, const Point &, update_parameters const& params) const {

void compute_reflection(Point &v, Point const&, update_parameters const& params) const {
Point a((-2.0 * params.inner_vi_ak) * A.row(params.facet_prev));
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
NT* v_data = v.pointerToData();
NT* p_data = p.pointerToData();
for(Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, params.facet_prev); it; ++it) {
*(v_data + it.col()) += (-2.0 * params.inner_vi_ak) * it.value();
*(p_data + it.col()) -= (-2.0 * params.inner_vi_ak * params.moved_dist) * it.value();
}
}

template <typename update_parameters>
NT compute_reflection(Point &v, const Point &, DenseMT const &AE, VT const &AEA, NT const &vEv, update_parameters const &params) const {

Expand Down
20 changes: 10 additions & 10 deletions include/preprocess/rounding_util_functions.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ initialize_chol(MT const& mat)
if constexpr (std::is_same<MT, DenseMT>::value)
{
return std::make_unique<Eigen::LLT<MT>>();
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
auto llt = std::make_unique<Eigen::SimplicialLLT<MT>>();
llt->analyzePattern(mat);
Expand All @@ -78,7 +78,7 @@ initialize_chol(MT const& A_trans, MT const& A)
if constexpr (std::is_same<MT, DenseMT>::value)
{
return std::make_unique<Eigen::LLT<MT>>();
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
MT mat = A_trans * A;
return initialize_chol<NT>(mat);
Expand All @@ -99,7 +99,7 @@ inline static VT solve_vec(std::unique_ptr<Eigen_lltMT> const& llt,
{
llt->compute(H);
return llt->solve(b);
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
llt->factorize(H);
return llt->solve(b);
Expand All @@ -122,7 +122,7 @@ solve_mat(std::unique_ptr<Eigen_lltMT> const& llt,
llt->compute(H);
logdetE = llt->matrixL().toDenseMatrix().diagonal().array().log().sum();
return llt->solve(mat);
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
llt->factorize(H);
logdetE = llt->matrixL().nestedExpression().diagonal().array().log().sum();
Expand All @@ -143,7 +143,7 @@ inline static void update_Atrans_Diag_A(MT &H, MT const& A_trans,
if constexpr (std::is_same<MT, DenseMT>::value)
{
H.noalias() = A_trans * D * A;
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
H = A_trans * D * A;
} else
Expand All @@ -161,7 +161,7 @@ inline static void update_Diag_A(MT &H, diag_MT const& D, MT const& A)
if constexpr (std::is_same<MT, DenseMT>::value)
{
H.noalias() = D * A;
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
H = D * A;
} else
Expand All @@ -179,7 +179,7 @@ inline static void update_A_Diag(MT &H, MT const& A, diag_MT const& D)
if constexpr (std::is_same<MT, DenseMT>::value)
{
H.noalias() = A * D;
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
H = A * D;
} else
Expand All @@ -198,7 +198,7 @@ get_mat_prod_op(MT const& E)
if constexpr (std::is_same<MT, DenseMT>::value)
{
return std::make_unique<Spectra::DenseSymMatProd<NT>>(E);
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
return std::make_unique<Spectra::SparseSymMatProd<NT>>(E);
} else
Expand Down Expand Up @@ -249,7 +249,7 @@ init_Bmat(MT &B, int const n, MT const& A_trans, MT const& A)
if constexpr (std::is_same<MT, DenseMT>::value)
{
B.resize(n+1, n+1);
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
// Initialize the structure of matrix B
typedef Eigen::Triplet<NT> triplet;
Expand Down Expand Up @@ -299,7 +299,7 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d,
B.block(n, 0, 1, n).noalias() = AtDe.transpose();
B(n, n) = d.sum();
B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1);
} else if constexpr (std::is_same<MT, SparseMT>::value)
} else if constexpr (std::is_base_of<Eigen::SparseMatrixBase<MT>, MT >::value)
{
MT AtD_A = AtD * A;
int k = 0;
Expand Down
2 changes: 1 addition & 1 deletion include/random_walks/gaussian_accelerated_billiard_walk.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,7 +73,7 @@ struct GaussianAcceleratedBilliardWalk

void computeLcov(const E_type E)
{
if constexpr (std::is_same<E_type, Eigen::SparseMatrix<NT> >::value) {
if constexpr (std::is_base_of<Eigen::SparseMatrixBase<E_type>, E_type >::value) {
Eigen::SimplicialLLT<E_type> lltofE;
lltofE.compute(E);
if (lltofE.info() != Eigen::Success) {
Expand Down
Loading
Loading