Skip to content

Commit

Permalink
New interface for sampling. Function is parameterized by walk and dis…
Browse files Browse the repository at this point in the history
…tribution. Added support for basic uniform samplers.
  • Loading branch information
vfisikop committed Apr 12, 2024
1 parent 6a5a17e commit 2527d7d
Show file tree
Hide file tree
Showing 7 changed files with 321 additions and 19 deletions.
29 changes: 29 additions & 0 deletions examples/sample_points/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
141 changes: 141 additions & 0 deletions examples/sample_points/simple_example.cpp
Original file line number Diff line number Diff line change
@@ -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<double> Kernel;
typedef typename Kernel::Point Point;
typedef BoostRandomNumberGenerator<boost::mt19937, double> RNGType;
typedef HPolytope<Point> HPolytopeType;

template <typename Walk, typename Distribution>
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<NT,Eigen::Dynamic,Eigen::Dynamic>;
using VT = Eigen::Matrix<NT,Eigen::Dynamic,1>;

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<NT, VT, MT>(samples, min_ess);
std::cout << "ess=" << min_ess << std::endl;
//print_diagnostics<NT, VT, MT>(samples, min_ess, std::cerr);
}

int main() {
// Generating a 3-dimensional cube centered at origin
HPolytopeType HP = generate_cube<HPolytopeType>(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<NT,Eigen::Dynamic,Eigen::Dynamic>;
using VT = Eigen::Matrix<NT,Eigen::Dynamic,1>;

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<Point> 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<HPolytopeType, RNGType> BilliardWalkType;
typedef AcceleratedBilliardWalk::template Walk<HPolytopeType, RNGType> AcceleratedBilliardWalkType;
typedef RandomPointGenerator<AcceleratedBilliardWalkType> Generator;
std::vector<Point> 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<Point> randPoints;
// gaussian_sampling<AcceleratedBilliardWalk>(points, HP, rng, walkL, numpoints, 1.0,
// StartingPoint, nburns);
double variance = 1.0;
Point c(HP.dimension());
HP.set_InnerBall(std::pair<Point,double>(Point(HP.dimension()), 1.0));
c = GetDirection<Point>::apply(HP.dimension(), rng, false);
//ExponentialHamiltonianMonteCarloExactWalk
exponential_sampling<NutsHamiltonianMonteCarloWalk>(points, HP, rng, walkL, numpoints, c, variance,
StartingPoint, nburns);
*/
return 0;
}
14 changes: 7 additions & 7 deletions include/random_walks/uniform_accelerated_billiard_walk.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -177,15 +177,15 @@ 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<Point>::apply(P.dimension(), radius, rng);
p += center;
initialize(P, p, rng);

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;
}
Expand All @@ -196,7 +196,7 @@ struct AcceleratedBilliardWalk
}

if (Lmax > _L) {
if (P.dimension() <= 2500)
if (P.dimension() <= 2500)
{
update_delta(Lmax);
}
Expand All @@ -207,7 +207,7 @@ struct AcceleratedBilliardWalk
pointset.clear();
}



inline void update_delta(NT L)
{
Expand Down Expand Up @@ -273,11 +273,11 @@ struct AcceleratedBilliardWalk
}
}

inline double get_max_distance(std::vector<Point> &pointset, Point const& q, double &rad)
inline double get_max_distance(std::vector<Point> &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) {
Expand Down
10 changes: 6 additions & 4 deletions include/random_walks/uniform_dikin_walk.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,20 +57,22 @@ struct DikinWalk
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;

Walk(Polytope &P, Point &p, RandomNumberGenerator &)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Expand Down
10 changes: 6 additions & 4 deletions include/random_walks/uniform_john_walk.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,20 +49,22 @@ struct JohnWalk
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;

Walk(Polytope &P, Point &p, RandomNumberGenerator &)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Expand Down
10 changes: 6 additions & 4 deletions include/random_walks/uniform_vaidya_walk.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,20 +50,22 @@ struct VaidyaWalk
typedef typename Polytope::VT VT;
typedef typename Point::FT NT;

Walk(Polytope &P, Point &p, RandomNumberGenerator &)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Walk(Polytope &P, Point & p, RandomNumberGenerator &, parameters const& params)
template <typename GenericPolytope>
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<NT>(p0, A, b, r);
}

Expand Down
Loading

0 comments on commit 2527d7d

Please sign in to comment.