From add54f40ca0d29a5eb21d52d254b201513b0ed4a Mon Sep 17 00:00:00 2001 From: vfisikop Date: Fri, 12 Apr 2024 16:35:37 +0300 Subject: [PATCH] New interface for sampling. Function is parameterized by walk and distribution. Added support for basic uniform samplers. --- examples/sample_points/CMakeLists.txt | 29 ++++ examples/sample_points/simple_example.cpp | 141 ++++++++++++++++++ .../uniform_accelerated_billiard_walk.hpp | 14 +- include/random_walks/uniform_dikin_walk.hpp | 10 +- include/random_walks/uniform_john_walk.hpp | 10 +- include/random_walks/uniform_vaidya_walk.hpp | 10 +- include/sampling/sample_points.hpp | 126 ++++++++++++++++ 7 files changed, 321 insertions(+), 19 deletions(-) create mode 100644 examples/sample_points/CMakeLists.txt create mode 100644 examples/sample_points/simple_example.cpp create mode 100644 include/sampling/sample_points.hpp diff --git a/examples/sample_points/CMakeLists.txt b/examples/sample_points/CMakeLists.txt new file mode 100644 index 000000000..bffa63a84 --- /dev/null +++ b/examples/sample_points/CMakeLists.txt @@ -0,0 +1,29 @@ +project( VolEsti-cpp-example ) + +CMAKE_MINIMUM_REQUIRED(VERSION 3.11) + +add_definitions(-DDISABLE_NLP_ORACLES) + +include("../../external/cmake-files/Eigen.cmake") +GetEigen() +if (${CMAKE_VERSION} VERSION_LESS "3.12.0") + add_compile_options(-D "EIGEN_NO_DEBUG") +else () + add_compile_definitions("EIGEN_NO_DEBUG") +endif () + +include("../../external/cmake-files/Boost.cmake") +GetBoost() + +include("../../external/cmake-files/LPSolve.cmake") +GetLPSolve() + +set(CMAKE_EXPORT_COMPILE_COMMANDS "ON") + +add_definitions(${CMAKE_CXX_FLAGS} "-g") + +include_directories (BEFORE ../../external) +include_directories (BEFORE ../../include) + +add_executable (example simple_example.cpp) +target_link_libraries(example PUBLIC lp_solve) \ No newline at end of file diff --git a/examples/sample_points/simple_example.cpp b/examples/sample_points/simple_example.cpp new file mode 100644 index 000000000..41fb28290 --- /dev/null +++ b/examples/sample_points/simple_example.cpp @@ -0,0 +1,141 @@ +#include "Eigen/Eigen" +#include "cartesian_geom/cartesian_kernel.h" +#include "convex_bodies/hpolytope.h" + +#include "diagnostics/diagnostics.hpp" +//todo make headers automomous +//#include "diagnostics/effective_sample_size.hpp" +//#include "diagnostics/print_diagnostics.hpp" + +#include "generators/known_polytope_generators.h" + +#include "random_walks/random_walks.hpp" + +#include "sampling/sample_points.hpp" +#include "sampling/sampling.hpp" + +#include "volume/volume_cooling_balls.hpp" + +typedef Cartesian Kernel; +typedef typename Kernel::Point Point; +typedef BoostRandomNumberGenerator RNGType; +typedef HPolytope HPolytopeType; + +template +void sample_points_eigen_matrix(HPolytopeType const& HP, Point const& q, Walk const& walk, + Distribution const& distr, RNGType rng, int walk_len, int rnum, + int nburns) +{ + using NT = double; + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + + MT samples(HP.dimension(), rnum); + + sample_points(HP, q, walk, distr, rng, walk_len, rnum, nburns, samples); + + // sample stats + unsigned int min_ess; + auto score = effective_sample_size(samples, min_ess); + std::cout << "ess=" << min_ess << std::endl; + //print_diagnostics(samples, min_ess, std::cerr); +} + +int main() { + // Generating a 3-dimensional cube centered at origin + HPolytopeType HP = generate_cube(10, false); + std::cout<<"Polytope: \n"; + HP.ComputeInnerBall(); + //HP.print(); + //std::cout<<"\n"; + + // Setup parameters for sampling + Point q(HP.dimension()); + RNGType rng(HP.dimension()); + + // NEW INTERFACE Sampling + AcceleratedBilliardWalk abill_walk; + AcceleratedBilliardWalk abill_walk_custom(10); //user defined walk parameters + BallWalk ball_walk; + BilliardWalk bill_walk; + CDHRWalk cdhr_walk; + DikinWalk dikin_walk; + JohnWalk john_walk; + RDHRWalk rdhr_walk; + VaidyaWalk vaidya_walk; + + GaussianBallWalk gball_walk; + + UniformDistribution udistr{}; + GaussianDistribution gdistr{}; + + using NT = double; + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + + int rnum = 100; + int nburns = 50; + int walk_len = 10; + + MT samples(HP.dimension(), rnum); + + // 1. the eigen matrix interface + sample_points_eigen_matrix(HP, q, abill_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, abill_walk_custom, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, ball_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, cdhr_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, dikin_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, john_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, rdhr_walk, udistr, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, vaidya_walk, udistr, rng, walk_len, rnum, nburns); + + // TODO: fix + // Does not converge because of the starting point + // Also ess returns rnum instead of 0 + sample_points_eigen_matrix(HP, q, bill_walk, udistr, rng, walk_len, rnum, nburns); + + // Does not compile because of walk-distribution combination + //sample_points_eigen_matrix(HP, q, abill_walk, gdistr, rng, walk_len, rnum, nburns); + + // 2. the std::vector interface + std::vector points; + sample_points(HP, q, cdhr_walk, udistr, rng, walk_len, rnum, nburns, points); + for (auto& point : points) + { + // std::cout << point.getCoefficients().transpose() << "\n"; + } + + // 3. the old interface + // different billiard walks + typedef BilliardWalk::template Walk BilliardWalkType; + typedef AcceleratedBilliardWalk::template Walk AcceleratedBilliardWalkType; + typedef RandomPointGenerator Generator; + std::vector randPoints; + PushBackWalkPolicy push_back_policy; + Generator::apply(HP, q, rnum, walk_len, randPoints, push_back_policy, rng); + for (auto& point : randPoints) + { + // std::cout << point.getCoefficients().transpose() << "\n"; + } + + +/* + unsigned int walkL = 10, numpoints = 10000, nburns = 0, d = HP.dimension(); + Point StartingPoint(d); + std::list randPoints; + +// gaussian_sampling(points, HP, rng, walkL, numpoints, 1.0, + // StartingPoint, nburns); + + double variance = 1.0; + + Point c(HP.dimension()); + HP.set_InnerBall(std::pair(Point(HP.dimension()), 1.0)); + c = GetDirection::apply(HP.dimension(), rng, false); + //ExponentialHamiltonianMonteCarloExactWalk + exponential_sampling(points, HP, rng, walkL, numpoints, c, variance, + StartingPoint, nburns); + +*/ + return 0; +} \ No newline at end of file diff --git a/include/random_walks/uniform_accelerated_billiard_walk.hpp b/include/random_walks/uniform_accelerated_billiard_walk.hpp index 980c8e5da..2c727a15c 100644 --- a/include/random_walks/uniform_accelerated_billiard_walk.hpp +++ b/include/random_walks/uniform_accelerated_billiard_walk.hpp @@ -165,7 +165,7 @@ struct AcceleratedBilliardWalk < typename GenericPolytope > - inline void parameters_burnin(GenericPolytope &P, + inline void parameters_burnin(GenericPolytope &P, Point const& center, unsigned int const& num_points, unsigned int const& walk_length, @@ -177,7 +177,7 @@ struct AcceleratedBilliardWalk pointset.push_back(_p); NT rad = NT(0), max_dist, Lmax = get_delta(), radius = P.InnerBall().second; - for (int i = 0; i < num_points; i++) + for (int i = 0; i < num_points; i++) { Point p = GetPointInDsphere::apply(P.dimension(), radius, rng); p += center; @@ -185,7 +185,7 @@ struct AcceleratedBilliardWalk apply(P, p, walk_length, rng); max_dist = get_max_distance(pointset, p, rad); - if (max_dist > Lmax) + if (max_dist > Lmax) { Lmax = max_dist; } @@ -196,7 +196,7 @@ struct AcceleratedBilliardWalk } if (Lmax > _L) { - if (P.dimension() <= 2500) + if (P.dimension() <= 2500) { update_delta(Lmax); } @@ -207,7 +207,7 @@ struct AcceleratedBilliardWalk pointset.clear(); } - + inline void update_delta(NT L) { @@ -273,11 +273,11 @@ struct AcceleratedBilliardWalk } } - inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) + inline double get_max_distance(std::vector &pointset, Point const& q, double &rad) { double dis = -1.0, temp_dis; int jj = 0; - for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) + for (auto vecit = pointset.begin(); vecit!=pointset.end(); vecit++, jj++) { temp_dis = (q.getCoefficients() - (*vecit).getCoefficients()).norm(); if (temp_dis > dis) { diff --git a/include/random_walks/uniform_dikin_walk.hpp b/include/random_walks/uniform_dikin_walk.hpp index e791b2544..326413d2c 100644 --- a/include/random_walks/uniform_dikin_walk.hpp +++ b/include/random_walks/uniform_dikin_walk.hpp @@ -57,20 +57,22 @@ struct DikinWalk typedef typename Polytope::VT VT; typedef typename Point::FT NT; - Walk(Polytope &P, Point &p, RandomNumberGenerator &) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); - NT r = P.ComputeInnerBall().second; + NT r = P.InnerBall().second; dikinw = DikinWalker(p0, A, b, r); } - Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); NT r = params.set_L ? params.m_L - : P.ComputeInnerBall().second; + : P.InnerBall().second; dikinw = DikinWalker(p0, A, b, r); } diff --git a/include/random_walks/uniform_john_walk.hpp b/include/random_walks/uniform_john_walk.hpp index 2c6bc42e3..cde92587a 100644 --- a/include/random_walks/uniform_john_walk.hpp +++ b/include/random_walks/uniform_john_walk.hpp @@ -49,20 +49,22 @@ struct JohnWalk typedef typename Polytope::VT VT; typedef typename Point::FT NT; - Walk(Polytope &P, Point &p, RandomNumberGenerator &) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); - NT r = P.ComputeInnerBall().second; + NT r = P.InnerBall().second; johnw = JohnWalker(p0, A, b, r); } - Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); NT r = params.set_L ? params.m_L - : P.ComputeInnerBall().second; + : P.InnerBall().second; johnw = JohnWalker(p0, A, b, r); } diff --git a/include/random_walks/uniform_vaidya_walk.hpp b/include/random_walks/uniform_vaidya_walk.hpp index 5bbbe9912..fd076cc28 100644 --- a/include/random_walks/uniform_vaidya_walk.hpp +++ b/include/random_walks/uniform_vaidya_walk.hpp @@ -50,20 +50,22 @@ struct VaidyaWalk typedef typename Polytope::VT VT; typedef typename Point::FT NT; - Walk(Polytope &P, Point &p, RandomNumberGenerator &) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); - NT r = P.ComputeInnerBall().second; + NT r = P.InnerBall().second; vaidyaw = VaidyaWalker(p0, A, b, r); } - Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params) + template + Walk(GenericPolytope& P, Point const& p, RandomNumberGenerator &, parameters const& params) { MT A = P.get_mat(); VT b = P.get_vec(), _vec_point = VT::Zero(P.dimension()), p0 = p.getCoefficients(); NT r = params.set_L ? params.m_L - : P.ComputeInnerBall().second; + : P.InnerBall().second; vaidyaw = VaidyaWalker(p0, A, b, r); } diff --git a/include/sampling/sample_points.hpp b/include/sampling/sample_points.hpp new file mode 100644 index 000000000..551316020 --- /dev/null +++ b/include/sampling/sample_points.hpp @@ -0,0 +1,126 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2024 Vissarion Fisikopoulos + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef SAMPLE_POINTS_HPP +#define SAMPLE_POINTS_HPP + +struct UniformDistribution {}; +struct GaussianDistribution {}; + +namespace detail +{ + +template +struct AlwaysFalse : std::false_type {}; + +template +struct sample_points +{ + static_assert(AlwaysFalse::value, "Not implemented for this combination of walk and distribution."); +}; + +} + +template +< + typename Polytope, + typename Point, + typename WalkType, + typename Distribution, + typename RandomNumberGenerator, + typename PointList +> +void sample_points(Polytope const& P, + Point const& starting_point, + WalkType const& walk_with_parameters, + Distribution const& distribution, + RandomNumberGenerator& rng, + unsigned int const& walk_len, + unsigned int const& rnum, + unsigned int const& nburns, + PointList& samples) +{ + if constexpr ((std::is_same::value + || std::is_same::value + || std::is_same::value + || std::is_same::value + || std::is_same::value + || std::is_same::value + || std::is_same::value + || std::is_same::value) + && std::is_same::value) + { + typename WalkType::template Walk + walk(P, starting_point, rng, walk_with_parameters.param); + + Point p = starting_point; + + // accelerated billiard walk variants in some cases need to compute a good starting point + // otherwise the walk stuck in the first given point (e.g. with a cube and its center of mass) + if constexpr (std::is_same::value) + walk.template get_starting_point(P, p, p, 10, rng); + + for (unsigned int i = 0; i < nburns; ++i) + { + walk.apply(P, p, walk_len, rng); + } + + for (unsigned int i = 0; i < rnum; ++i) + { + walk.apply(P, p, walk_len, rng); + samples.push_back(p); + } + } + else if constexpr (std::is_same::value + && std::is_same::value) + { + //TODO + } + else + { + static_assert(detail::AlwaysFalse::value, + "Not implemented for this combination of walk and distribution."); + } +} + + +// this is a wrapper function for storing the samples in an Eigen matrix +// it is inefficient in many ways (e.g. unneeded copies) and should be rewritten +template +< + typename Polytope, + typename Point, + typename WalkType, + typename Distribution, + typename RandomNumberGenerator, + typename NT +> +void sample_points(Polytope const& P, + Point const& starting_point, + WalkType const& walk_with_parameters, + Distribution const& distribution, + RandomNumberGenerator& rng, + unsigned int const& walk_len, + unsigned int const& rnum, + unsigned int const& nburns, + Eigen::Matrix& samples) +{ + std::vector random_point(1); + auto p = starting_point; + sample_points(P, p, walk_with_parameters, distribution, rng, walk_len, 0, nburns, random_point); + random_point.clear(); + + for (unsigned int i = 0; i < rnum; ++i) + { + sample_points(P, p, walk_with_parameters, distribution, rng, walk_len, 1, 0, random_point); + auto point = random_point[0]; + samples.col(i) = point.getCoefficients(); + random_point.clear(); + } + +} + +#endif //SAMPLE_POINTS_HPP \ No newline at end of file