From 8dc74d7d88733b4606a71aef1927449cce6915f9 Mon Sep 17 00:00:00 2001 From: vgnecula Date: Mon, 8 Apr 2024 02:42:37 -0400 Subject: [PATCH] Efficient BurnIn Updated --- include/sampling/sampling.hpp | 657 +++++++++++++++------------------- 1 file changed, 297 insertions(+), 360 deletions(-) diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp index 696186098..4d158f530 100644 --- a/include/sampling/sampling.hpp +++ b/include/sampling/sampling.hpp @@ -22,45 +22,41 @@ #ifndef SAMPLE_ONLY_H #define SAMPLE_ONLY_H -template +template +< + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename Point +> void uniform_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + Point &starting_point, + unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + typename WalkTypePolicy::template Walk walk(P, starting_point, rng); Point p = starting_point; - typedef RandomPointGenerator RandomPointGenerator; if (nburns > 0) { - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walk.apply(P, p, walk_len, rng); + } } - RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, - push_back_policy, rng); - + for (unsigned int i = 0; i < rnum; ++i) { + walk.apply(P, p, walk_len, rng); + randPoints.push_back(p); + } } -template < +template +< typename PointList, typename Polytope, typename RandomNumberGenerator, @@ -68,83 +64,76 @@ template < typename Point > void uniform_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - WalkTypePolicy &WalkType, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + WalkTypePolicy &WalkType, + const unsigned int &walk_len, + const unsigned int &rnum, + Point &starting_point, + unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; - typedef RandomPointGenerator RandomPointGenerator; + typename WalkTypePolicy::template Walk walk(P, starting_point, rng); Point p = starting_point; + if (nburns > 0) { - RandomPointGenerator::apply(P, p, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + 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); + randPoints.push_back(p); } - RandomPointGenerator::apply(P, p, rnum, walk_len, randPoints, - push_back_policy, rng, WalkType.param); } -template +template < - typename WalkTypePolicy, - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename Point + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename Point > void uniform_sampling_boundary(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &starting_point, + unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; + typedef typename WalkTypePolicy::template Walk Walk; - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; - - Point p = starting_point; + Walk walk(P, starting_point, rng); + Point current_point = starting_point; + Point dummy_point(P.dimension()); - typedef BoundaryRandomPointGenerator BoundaryRandomPointGenerator; - if (nburns > 0) { - BoundaryRandomPointGenerator::apply(P, p, nburns, walk_len, - randPoints, push_back_policy, rng); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walk.apply(P, current_point, dummy_point, walk_len, rng); + current_point = dummy_point; } - unsigned int n = rnum / 2; - BoundaryRandomPointGenerator::apply(P, p, rnum / 2, walk_len, - randPoints, push_back_policy, rng); + for (unsigned int i = 0; i < rnum; ++i) { + walk.apply(P, current_point, dummy_point, walk_len, rng); + randPoints.push_back(current_point); + current_point = dummy_point; + } } -template + +template < - typename WalkTypePolicy, - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename NT, - typename Point + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename NT, + typename Point > void gaussian_sampling(PointList &randPoints, Polytope &P, @@ -152,85 +141,74 @@ void gaussian_sampling(PointList &randPoints, const unsigned int &walk_len, const unsigned int &rnum, const NT &a, - const Point &starting_point, - unsigned int const& nburns) + const Point &starting_point, + unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + typename WalkTypePolicy::template Walk walk(P, starting_point, a, rng); Point p = starting_point; - typedef GaussianRandomPointGenerator RandomPointGenerator; if (nburns > 0) { - RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walk.apply(P, p, a, walk_len, rng); + } } - RandomPointGenerator::apply(P, p, a, rnum, walk_len, randPoints, - push_back_policy, rng); - + for (unsigned int i = 0; i < rnum; ++i) { + walk.apply(P, p, a, walk_len, rng); + randPoints.push_back(p); + } } -template < - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename WalkTypePolicy, - typename NT, - typename Point - > +template +< + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point +> void gaussian_sampling(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, - WalkTypePolicy &WalkType, + WalkTypePolicy &WalkType, const unsigned int &walk_len, const unsigned int &rnum, const NT &a, const Point &starting_point, - unsigned int const& nburns) + unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + typename WalkTypePolicy::template Walk walk(P, starting_point, a, rng); Point p = starting_point; - - typedef GaussianRandomPointGenerator RandomPointGenerator; + if (nburns > 0) { - RandomPointGenerator::apply(P, p, a, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walk.apply(P, p, a, walk_len, rng); + } + } + + for (unsigned int i = 0; i < rnum; ++i) { + walk.apply(P, p, a, walk_len, rng); + randPoints.push_back(p); } - RandomPointGenerator::apply(P, p, a, rnum, walk_len, randPoints, - push_back_policy, rng, WalkType.param); } template < - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename WalkTypePolicy, - typename NT, - typename Point, - typename NegativeGradientFunctor, - typename NegativeLogprobFunctor, - typename Solver - > + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point, + typename NegativeGradientFunctor, + typename NegativeLogprobFunctor, + typename Solver +> void logconcave_sampling(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, @@ -239,65 +217,47 @@ void logconcave_sampling(PointList &randPoints, const Point &starting_point, unsigned int const& nburns, NegativeGradientFunctor &F, - NegativeLogprobFunctor &f) + NegativeLogprobFunctor &f) { - typedef typename WalkTypePolicy::template Walk - < - Point, - Polytope, - RandomNumberGenerator, - NegativeGradientFunctor, - NegativeLogprobFunctor, - Solver - > walk; - - typedef typename WalkTypePolicy::template parameters - < - NT, - NegativeGradientFunctor - > walk_params; - - // Initialize random walk parameters - unsigned int dim = starting_point.dimension(); - walk_params params(F, dim); - - if (F.params.eta > 0) { - params.eta = F.params.eta; - } + + typename WalkTypePolicy::template Walk< + Polytope, RandomNumberGenerator, NegativeGradientFunctor, NegativeLogprobFunctor, Solver + > logconcave_walk(P, starting_point, F, f, rng); - PushBackWalkPolicy push_back_policy; + Point current_point = starting_point; + for (unsigned int i = 0; i < nburns; ++i) { + logconcave_walk.apply(current_point, walk_len, rng); + } - Point p = starting_point; + typedef LogconcaveRandomPointGenerator> RandomPointGenerator; - walk logconcave_walk(&P, p, F, f, params); + PushBackWalkPolicy push_back_policy; - typedef LogconcaveRandomPointGenerator RandomPointGenerator; - - if (nburns > 0) { - RandomPointGenerator::apply(nburns, walk_len, randPoints, - push_back_policy, rng, logconcave_walk); + for (unsigned int i = 0; i < rnum; ++i) { + logconcave_walk.apply(current_point, walk_len, rng); + push_back_policy.apply(randPoints, current_point); } - logconcave_walk.disable_adaptive(); - randPoints.clear(); - - RandomPointGenerator::apply(rnum, walk_len, randPoints, - push_back_policy, rng, logconcave_walk); } + #include "preprocess/crhmc/crhmc_input.h" #include "preprocess/crhmc/crhmc_problem.h" -template - < - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename WalkTypePolicy, - typename NT, - typename Point, - typename NegativeGradientFunctor, - typename NegativeLogprobFunctor, - typename HessianFunctor, - typename Solver - > +#include "random_walks/crhmc/crhmc_walk.hpp" +#include "random_walks/crhmc/crhmc_walk.hpp" + +template < + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point, + typename NegativeGradientFunctor, + typename NegativeLogprobFunctor, + typename HessianFunctor, + typename Solver +> void crhmc_sampling(PointList &randPoints, Polytope &P, RandomNumberGenerator &rng, @@ -308,55 +268,51 @@ void crhmc_sampling(PointList &randPoints, NegativeLogprobFunctor &f, HessianFunctor &h, int simdLen = 1, - bool raw_output=false) { - typedef typename Polytope::MT MatrixType; - typedef crhmc_input - < - MatrixType, - Point, - NegativeLogprobFunctor, - NegativeGradientFunctor, - HessianFunctor - > Input; - Input input = convert2crhmc_input(P, f, F, h); - typedef crhmc_problem CrhmcProblem; - CrhmcProblem problem = CrhmcProblem(input); - if(problem.terminate){return;} - typedef typename WalkTypePolicy::template Walk - < - Point, - CrhmcProblem, - RandomNumberGenerator, - NegativeGradientFunctor, - NegativeLogprobFunctor, - Solver - > walk; - typedef typename WalkTypePolicy::template parameters - < - NT, - NegativeGradientFunctor - > walk_params; - Point p = Point(problem.center); - problem.options.simdLen=simdLen; - walk_params params(input.df, p.dimension(), problem.options); - - if (input.df.params.eta > 0) { - params.eta = input.df.params.eta; - } - - PushBackWalkPolicy push_back_policy; - - walk crhmc_walk = walk(problem, p, input.df, input.f, params); - - typedef CrhmcRandomPointGenerator RandomPointGenerator; - - RandomPointGenerator::apply(problem, p, nburns, walk_len, randPoints, - push_back_policy, rng, F, f, params, crhmc_walk); - //crhmc_walk.disable_adaptive(); - randPoints.clear(); - RandomPointGenerator::apply(problem, p, rnum, walk_len, randPoints, - push_back_policy, rng, F, f, params, crhmc_walk, simdLen, raw_output); + bool raw_output=false) +{ + + typedef typename Polytope::MT MatrixType; + typedef crhmc_input Input; + + Input input = convert2crhmc_input(P, f, F, h); + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + if(problem.terminate) { return; } + + typename CRHMCWalk::template parameters params(F, P.dimension(), problem.options); + + if (input.df.params.eta > 0) { + params.eta = input.df.params.eta; + } + + Point p = Point(problem.center); + + typename CRHMCWalk::template Walk< + Point, + CrhmcProblem, + RandomNumberGenerator, + NegativeGradientFunctor, + NegativeLogprobFunctor, + Solver + > crhmc_walk(problem, p, F, f, params); + + for(unsigned int i = 0; i < nburns; ++i) { + crhmc_walk.apply(rng, walk_len); + } + + for(unsigned int i = 0; i < rnum; ++i) { + crhmc_walk.apply(rng, walk_len); + if (raw_output) { + + } else { + randPoints.push_back(crhmc_walk.getPoint()); + } + } } + + + + #include "ode_solvers/ode_solvers.hpp" template < typename Polytope, @@ -434,169 +390,150 @@ crhmc_sampling < >(randPoints, P, rng, walkL, numpoints, nburns, *F, *f, zerof, simdLen, raw_output); } } -template -< - typename WalkTypePolicy, - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename NT, - typename Point + +template < + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename NT, + typename Point > void exponential_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &c, - const NT &a, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + Point &starting_point, + unsigned int const& nburns) { - - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; + typedef typename WalkTypePolicy::template Walk walk; PushBackWalkPolicy push_back_policy; + Point p = starting_point; - Point p = starting_point; + typedef ExponentialRandomPointGenerator RandomPointGenerator; - typedef ExponentialRandomPointGenerator RandomPointGenerator; if (nburns > 0) { - RandomPointGenerator::apply(P, p, c, a, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walk tmpWalk(P, p, rng); + tmpWalk.apply(p, c, a, walk_len); + } + } + + for (unsigned int i = 0; i < rnum; ++i) { + walk tmpWalk(P, p, rng); + tmpWalk.apply(p, c, a, walk_len); + randPoints.push_back(p); } - RandomPointGenerator::apply(P, p, c, a, rnum, walk_len, randPoints, - push_back_policy, rng); } template < - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename WalkTypePolicy, - typename NT, - typename Point - > + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point +> void exponential_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - WalkTypePolicy &WalkType, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &c, - const NT &a, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + WalkTypePolicy &WalkType, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + Point &starting_point, + unsigned int const& nburns) { + typedef typename WalkTypePolicy::template Walk walk; + + walk walkInstance(P, starting_point, rng, WalkType.param); - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - PushBackWalkPolicy push_back_policy; - - Point p = starting_point; - - typedef ExponentialRandomPointGenerator RandomPointGenerator; if (nburns > 0) { - RandomPointGenerator::apply(P, p, c, a, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walkInstance.apply(starting_point, c, a, walk_len); + } + } + for (unsigned int i = 0; i < rnum; ++i) { + walkInstance.apply(starting_point, c, a, walk_len); + randPoints.push_back(starting_point); } - RandomPointGenerator::apply(P, p, c, a, rnum, walk_len, randPoints, - push_back_policy, rng, WalkType.param); } - -template -< - typename WalkTypePolicy, - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename NT, - typename Point +template < + typename WalkTypePolicy, + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename NT, + typename Point > void exponential_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &c, - const NT &a, - const NT &eta, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const NT &eta, + Point &starting_point, + unsigned int const& nburns) { + typedef typename WalkTypePolicy::template Walk walk; - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - PushBackWalkPolicy push_back_policy; - - Point p = starting_point; + walk walkInstance(P, starting_point, rng); - typedef ExponentialRandomPointGenerator RandomPointGenerator; if (nburns > 0) { - RandomPointGenerator::apply(P, p, c, a, eta, nburns, walk_len, randPoints, - push_back_policy, rng); - randPoints.clear(); + for (unsigned int i = 0; i < nburns; ++i) { + walkInstance.apply(starting_point, c, a, eta, walk_len); + } + } +\ + for (unsigned int i = 0; i < rnum; ++i) { + walkInstance.apply(starting_point, c, a, eta, walk_len); + randPoints.push_back(starting_point); } - RandomPointGenerator::apply(P, p, c, a, eta, rnum, walk_len, randPoints, - push_back_policy, rng); } template < - typename PointList, - typename Polytope, - typename RandomNumberGenerator, - typename WalkTypePolicy, - typename NT, - typename Point - > + typename PointList, + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy, + typename NT, + typename Point +> void exponential_sampling(PointList &randPoints, - Polytope &P, - RandomNumberGenerator &rng, - WalkTypePolicy &WalkType, - const unsigned int &walk_len, - const unsigned int &rnum, - const Point &c, - const NT &a, - const NT &eta, - const Point &starting_point, - unsigned int const& nburns) + Polytope &P, + RandomNumberGenerator &rng, + WalkTypePolicy &WalkType, + const unsigned int &walk_len, + const unsigned int &rnum, + const Point &c, + const NT &a, + const NT &eta, + const Point &starting_point, + unsigned int const& nburns) { - - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - PushBackWalkPolicy push_back_policy; + + typename WalkTypePolicy::template Walk walk(P, starting_point, rng, WalkType.param); Point p = starting_point; - typedef ExponentialRandomPointGenerator RandomPointGenerator; - if (nburns > 0) { - RandomPointGenerator::apply(P, p, c, a, eta, nburns, walk_len, randPoints, - push_back_policy, rng, WalkType.param); - randPoints.clear(); + for(unsigned int i = 0; i < nburns; ++i) { + walk.apply(p, walk_len, c, a, eta, rng); + } + + for(unsigned int i = 0; i < rnum; ++i) { + walk.apply(p, walk_len, c, a, eta, rng); + randPoints.push_back(p); } - RandomPointGenerator::apply(P, p, c, a, eta, rnum, walk_len, randPoints, - push_back_policy, rng, WalkType.param); }