From 99cbba620ba40f534f77999273f9b2a1cad5dbf2 Mon Sep 17 00:00:00 2001 From: vfisikop Date: Fri, 5 Jul 2024 12:55:08 +0300 Subject: [PATCH] Simplify sample_points and use NT template instead of double in distribution classes --- examples/sample_points/simple_example.cpp | 7 ++- include/sampling/sample_points.hpp | 70 +++++++++++------------ 2 files changed, 39 insertions(+), 38 deletions(-) diff --git a/examples/sample_points/simple_example.cpp b/examples/sample_points/simple_example.cpp index 30bf154c9..27b496cb4 100644 --- a/examples/sample_points/simple_example.cpp +++ b/examples/sample_points/simple_example.cpp @@ -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(end - start).count(); @@ -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 diff --git a/include/sampling/sample_points.hpp b/include/sampling/sample_points.hpp index 8eaf2d662..bc6d341cf 100644 --- a/include/sampling/sample_points.hpp +++ b/include/sampling/sample_points.hpp @@ -11,28 +11,31 @@ #include "preprocess/crhmc/crhmc_problem.h" #include "sampling/sampling.hpp" +template struct UniformDistribution { - using CartesianPoint = point>; + using CartesianPoint = point>; using Func = ZeroScalarFunctor; using Grad = ZeroFunctor; using Hess = ZeroFunctor; }; +template struct SphericalGaussianDistribution { SphericalGaussianDistribution() : variance(1.0) {} - SphericalGaussianDistribution(double _variance) + SphericalGaussianDistribution(NT _variance) : variance(_variance) {} - double variance; + NT variance; }; +template struct GaussianDistribution { - using CartesianEllipsoid = Ellipsoid>>; + using CartesianEllipsoid = Ellipsoid>>; GaussianDistribution() {} @@ -42,45 +45,46 @@ struct GaussianDistribution CartesianEllipsoid ellipsoid; }; +template struct ExponentialDistribution { - using CartesianPoint = point>; + using CartesianPoint = point>; ExponentialDistribution() {} - ExponentialDistribution(CartesianPoint _c, double _T) + ExponentialDistribution(CartesianPoint _c, NT _T) : c(_c) , T(_T) {} CartesianPoint c; - double T; + NT T; }; - +template struct LogConcaveDistribution { - using CartesianPoint = point>; + using CartesianPoint = point>; template - LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, double _L) + LogConcaveDistribution(GradientFunctor g, FunctionFunctor f, NT _L) : grad(g) , func(f) , L(_L) {} template - 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 const&, double const&)> grad; + std::function const&, NT const&)> grad; std::function grad_point; - std::function func; + std::function func; std::function hess; - double L; + NT L; }; namespace detail @@ -99,6 +103,7 @@ struct sample_points template < + int simdLen = 1, typename Polytope, typename Point, typename WalkType, @@ -124,7 +129,7 @@ void sample_points(Polytope& P, // TODO: make it a const& || 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); @@ -146,7 +151,7 @@ void sample_points(Polytope& P, // TODO: make it a const& || 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, distribution.variance, rng, walk_with_parameters.param); @@ -165,7 +170,7 @@ void sample_points(Polytope& P, // TODO: make it a const& } } else if constexpr (std::is_same::value - && std::is_same::value) + && std::is_same>::value) { typename WalkType::template Walk walk(P, starting_point, distribution.ellipsoid, rng, walk_with_parameters.param); @@ -184,7 +189,7 @@ void sample_points(Polytope& P, // TODO: make it a const& } } else if constexpr (std::is_same::value - && std::is_same::value) + && std::is_same>::value) { typename WalkType::template Walk walk(P, starting_point, distribution.c, distribution.T, rng, walk_with_parameters.param); @@ -205,19 +210,19 @@ void sample_points(Polytope& P, // TODO: make it a const& else if constexpr ((std::is_same::value || std::is_same::value || std::is_same::value) - && std::is_same::value) + && std::is_same>::value) { using HPolytope = typename std::remove_const::type; - using Solver = LeapfrogODESolver; + using Solver = LeapfrogODESolver; std::vector 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()); @@ -244,13 +249,11 @@ void sample_points(Polytope& P, // TODO: make it a const& } } else if constexpr ((std::is_same::value) - && std::is_same::value) + && std::is_same>::value) { using HPolytope = typename std::remove_const::type; HPolytope HP = P; //TODO: avoid the copy - constexpr int simdLen = 8; //TODO: input parameter - int dimension = HP.dimension(); auto G = distribution.grad_point; @@ -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(HP, F, G, H); using CrhmcProblem = crhmc_problem; CrhmcProblem problem = CrhmcProblem(input); @@ -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 @@ -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; // } @@ -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) { @@ -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); } } }