Skip to content

Commit

Permalink
Simplify sample_points and use NT template instead of double in distr…
Browse files Browse the repository at this point in the history 
…ibution classes
  • Loading branch information
@vfisikop
vfisikop committed Sep 10, 2024
1 parent 459879a commit 99cbba6
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 38 deletions.
7 changes: 5 additions & 2 deletions examples/sample_points/simple_example.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@ void sample_points_eigen_matrix(PolytopeOrProblem const& HP, Point const& q, Wal

auto start = std::chrono::steady_clock::now();

sample_points(HP, q, walk, distr, rng, walk_len, rnum, nburns, samples);
sample_points<64>(HP, q, walk, distr, rng, walk_len, rnum, nburns, samples);

auto end = std::chrono::steady_clock::now();
auto time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
Expand Down Expand Up @@ -314,10 +314,13 @@ int main() {
// The following will compile but segfauls since walk and distribution are not compatible
//sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns);
sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns);
sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_uniform, rng, walk_len, rnum, nburns);
sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_uniform, rng, walk_len, 100000, nburns);
sample_points_eigen_matrix(HP, q, crhmc_walk, logconcave_uniform, rng, walk_len, rnum, nburns);
//sample_points_eigen_matrix(problem, q, nhmc_walk, logconcave_uniform2, rng, walk_len, rnum, nburns);

using NT_fromMT = VT::Scalar;
NT_fromMT test;

std::cout << "fix the following" << std::endl;
// TODO: fix
// Does not converge because of the starting point
Expand Down
70 changes: 34 additions & 36 deletions include/sampling/sample_points.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,28 +11,31 @@
#include "preprocess/crhmc/crhmc_problem.h"
#include "sampling/sampling.hpp"

template <typename NT = double>
struct UniformDistribution
{
using CartesianPoint = point<Cartesian<double>>;
using CartesianPoint = point<Cartesian<NT>>;
using Func = ZeroScalarFunctor<CartesianPoint>;
using Grad = ZeroFunctor<CartesianPoint>;
using Hess = ZeroFunctor<CartesianPoint>;
};

template <typename NT = double>
struct SphericalGaussianDistribution
{
SphericalGaussianDistribution()
: variance(1.0)
{}
SphericalGaussianDistribution(double _variance)
SphericalGaussianDistribution(NT _variance)
: variance(_variance)
{}
double variance;
NT variance;
};

template <typename NT = double>
struct GaussianDistribution
{
using CartesianEllipsoid = Ellipsoid<point<Cartesian<double>>>;
using CartesianEllipsoid = Ellipsoid<point<Cartesian<NT>>>;

GaussianDistribution() {}

Expand All @@ -42,45 +45,46 @@ struct GaussianDistribution
CartesianEllipsoid ellipsoid;
};

template <typename NT = double>
struct ExponentialDistribution
{
using CartesianPoint = point<Cartesian<double>>;
using CartesianPoint = point<Cartesian<NT>>;

ExponentialDistribution() {}

ExponentialDistribution(CartesianPoint _c, double _T)
ExponentialDistribution(CartesianPoint _c, NT _T)
: c(_c)
, T(_T)
{}
CartesianPoint c;
double T;
NT T;
};


template <typename NT = double>
struct LogConcaveDistribution
{
using CartesianPoint = point<Cartesian<double>>;
using CartesianPoint = point<Cartesian<NT>>;

template <typename GradientFunctor, typename FunctionFunctor>
LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, double _L)
LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, NT _L)
: grad(g)
, func(f)
, L(_L)
{}

template <typename GradientFunctor, typename FunctionFunctor, typename HessianFunctor>
LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, HessianFunctor h, double _L)
LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, HessianFunctor h, NT _L)
: grad_point(g)
, func(f)
, hess(h)
, L(_L)
{}

std::function<CartesianPoint(unsigned int const&, std::vector<CartesianPoint> const&, double const&)> grad;
std::function<CartesianPoint(unsigned int const&, std::vector<CartesianPoint> const&, NT const&)> grad;
std::function<CartesianPoint(CartesianPoint const&)> grad_point;
std::function<double(CartesianPoint const&)> func;
std::function<NT(CartesianPoint const&)> func;
std::function<CartesianPoint(CartesianPoint const&)> hess;
double L;
NT L;
};

namespace detail
Expand All @@ -99,6 +103,7 @@ struct sample_points

template
<
int simdLen = 1,
typename Polytope,
typename Point,
typename WalkType,
Expand All @@ -124,7 +129,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
|| std::is_same<WalkType, JohnWalk>::value
|| std::is_same<WalkType, RDHRWalk>::value
|| std::is_same<WalkType, VaidyaWalk>::value)
&& std::is_same<Distribution, UniformDistribution>::value)
&& std::is_same<Distribution, UniformDistribution<NT>>::value)
{
typename WalkType::template Walk<Polytope, RandomNumberGenerator>
walk(P, starting_point, rng, walk_with_parameters.param);
Expand All @@ -146,7 +151,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
|| std::is_same<WalkType, GaussianCDHRWalk>::value
|| std::is_same<WalkType, GaussianHamiltonianMonteCarloExactWalk>::value
|| std::is_same<WalkType, GaussianRDHRWalk>::value)
&& std::is_same<Distribution, SphericalGaussianDistribution>::value)
&& std::is_same<Distribution, SphericalGaussianDistribution<NT>>::value)
{
typename WalkType::template Walk<Polytope, RandomNumberGenerator>
walk(P, starting_point, distribution.variance, rng, walk_with_parameters.param);
Expand All @@ -165,7 +170,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
}
}
else if constexpr (std::is_same<WalkType, GaussianAcceleratedBilliardWalk>::value
&& std::is_same<Distribution, GaussianDistribution>::value)
&& std::is_same<Distribution, GaussianDistribution<NT>>::value)
{
typename WalkType::template Walk<Polytope, RandomNumberGenerator>
walk(P, starting_point, distribution.ellipsoid, rng, walk_with_parameters.param);
Expand All @@ -184,7 +189,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
}
}
else if constexpr (std::is_same<WalkType, ExponentialHamiltonianMonteCarloExactWalk>::value
&& std::is_same<Distribution, ExponentialDistribution>::value)
&& std::is_same<Distribution, ExponentialDistribution<NT>>::value)
{
typename WalkType::template Walk<Polytope, RandomNumberGenerator>
walk(P, starting_point, distribution.c, distribution.T, rng, walk_with_parameters.param);
Expand All @@ -205,19 +210,19 @@ void sample_points(Polytope& P, // TODO: make it a const&
else if constexpr ((std::is_same<WalkType, HamiltonianMonteCarloWalk>::value
|| std::is_same<WalkType, NutsHamiltonianMonteCarloWalk>::value
|| std::is_same<WalkType, UnderdampedLangevinWalk>::value)
&& std::is_same<Distribution, LogConcaveDistribution>::value)
&& std::is_same<Distribution, LogConcaveDistribution<NT>>::value)
{
using HPolytope = typename std::remove_const<Polytope>::type;

using Solver = LeapfrogODESolver<Point, double, HPolytope, decltype(distribution.grad)>;
using Solver = LeapfrogODESolver<Point, NT, HPolytope, decltype(distribution.grad)>;

std::vector<Point> xs;
unsigned int i = 0;
double t = 1.0;
NT t = 1.0;

typename WalkType::template parameters
<
double,
NT,
decltype(distribution.grad)
> hmc_params(distribution.L, P.dimension());

Expand All @@ -244,13 +249,11 @@ void sample_points(Polytope& P, // TODO: make it a const&
}
}
else if constexpr ((std::is_same<WalkType, CRHMCWalk>::value)
&& std::is_same<Distribution, LogConcaveDistribution>::value)
&& std::is_same<Distribution, LogConcaveDistribution<NT>>::value)
{
using HPolytope = typename std::remove_const<Polytope>::type;
HPolytope HP = P; //TODO: avoid the copy

constexpr int simdLen = 8; //TODO: input parameter

int dimension = HP.dimension();

auto G = distribution.grad_point;
Expand All @@ -270,11 +273,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
NegativeGradientFunctor,
HessianFunctor
>;
Input input = convert2crhmc_input
<
Input, HPolytope, NegativeLogprobFunctor,
NegativeGradientFunctor, HessianFunctor
>(HP, F, G, H);
Input input = convert2crhmc_input<Input>(HP, F, G, H);

using CrhmcProblem = crhmc_problem<Point, Input>;
CrhmcProblem problem = CrhmcProblem(input);
Expand All @@ -283,7 +282,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
using Solver = ImplicitMidpointODESolver
<
Point, NT, CrhmcProblem,
decltype(distribution.grad_point), simdLen
NegativeGradientFunctor, simdLen
>;

using Walk = typename WalkType::template Walk
Expand All @@ -301,10 +300,10 @@ void sample_points(Polytope& P, // TODO: make it a const&
NegativeGradientFunctor
>;
Point p = Point(problem.center);
problem.options.simdLen=simdLen;
problem.options.simdLen = simdLen;
WalkParams params(distribution.L, p.dimension(), problem.options);

// TODO: pass a usef defined eta to bypass the one created by the walk using L
// TODO: pass a user defined eta to bypass the one created by the walk using L
//if (distribution.eta > 0) {
// params.eta = distribution.eta;
// }
Expand All @@ -316,12 +315,11 @@ void sample_points(Polytope& P, // TODO: make it a const&
walk.apply(rng, walk_len);
}

bool raw_output = false; // TODO: check this
for (unsigned int i = 0; i < std::ceil((float)rnum/simdLen); ++i)
{
walk.apply(rng, walk_len);
if (walk.P.terminate) {return;}
auto x = raw_output ? walk.x : walk.getPoints();
auto x = walk.getPoints();

if ((i + 1) * simdLen > rnum)
{
Expand All @@ -333,7 +331,7 @@ void sample_points(Polytope& P, // TODO: make it a const&
}
for (int j = 0; j < x.cols(); j++)
{
samples.col(i+j) = x.col(j);
samples.col(i+j) = x.col(j);
}
}
}
Expand Down

0 comments on commit 99cbba6

Please sign in to comment.