From 8367fa338df6d22d32df9d1e1e89eeb975ce54c6 Mon Sep 17 00:00:00 2001 From: vgnecula Date: Sat, 20 Apr 2024 01:30:30 -0400 Subject: [PATCH] BurnIn Handling, same style --- include/sampling/sampling.hpp | 351 ++++++++++++++++------------------ 1 file changed, 167 insertions(+), 184 deletions(-) diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp index 696186098..708e73df6 100644 --- a/include/sampling/sampling.hpp +++ b/include/sampling/sampling.hpp @@ -33,31 +33,28 @@ void uniform_sampling(PointList &randPoints, RandomNumberGenerator &rng, const unsigned int &walk_len, const unsigned int &rnum, - const Point &starting_point, + Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk + typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + > 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 < @@ -73,27 +70,28 @@ void uniform_sampling(PointList &randPoints, WalkTypePolicy &WalkType, const unsigned int &walk_len, const unsigned int &rnum, - const Point &starting_point, + Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk + + typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; - typedef RandomPointGenerator RandomPointGenerator; + > 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); } @@ -115,25 +113,24 @@ void uniform_sampling_boundary(PointList &randPoints, { typedef typename WalkTypePolicy::template Walk < - Polytope, - RandomNumberGenerator - > walk; + Polytope, + RandomNumberGenerator + > Walk; - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + Walk walk(P, starting_point, rng); + Point current_point = starting_point; + Point dummy_point(P.dimension()); - Point p = starting_point; - - 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/2; ++i) { + walk.apply(P, current_point, dummy_point, walk_len, rng); + randPoints.push_back(current_point); + current_point = dummy_point; + } } @@ -156,27 +153,24 @@ void gaussian_sampling(PointList &randPoints, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk + typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + > 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); + } } @@ -199,25 +193,24 @@ void gaussian_sampling(PointList &randPoints, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < + typename WalkTypePolicy::template Walk + < Polytope, RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + > 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 < @@ -241,50 +234,38 @@ void logconcave_sampling(PointList &randPoints, NegativeGradientFunctor &F, NegativeLogprobFunctor &f) { - typedef typename WalkTypePolicy::template Walk + + typename WalkTypePolicy::template Walk < - Point, Polytope, RandomNumberGenerator, NegativeGradientFunctor, NegativeLogprobFunctor, Solver - > walk; - - typedef typename WalkTypePolicy::template parameters - < - NT, - NegativeGradientFunctor - > walk_params; + > logconcave_walk(P, starting_point, F, f, rng); - // Initialize random walk parameters - unsigned int dim = starting_point.dimension(); - walk_params params(F, dim); + auto params = typename WalkTypePolicy::parameters(F, starting_point.dimension()); + logconcave_walk.set_parameters(params); - if (F.params.eta > 0) { - params.eta = F.params.eta; + Point current_point = starting_point; + for (unsigned int i = 0; i < nburns; ++i) { + logconcave_walk.apply(current_point, walk_len, rng); } PushBackWalkPolicy push_back_policy; - - Point p = starting_point; - - walk logconcave_walk(&P, p, F, f, params); - - 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" +#include "random_walks/crhmc/crhmc_walk.hpp" +#include "random_walks/crhmc/crhmc_walk.hpp" + template < typename PointList, @@ -309,54 +290,59 @@ void crhmc_sampling(PointList &randPoints, 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); + 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; } + + Point p = Point(problem.center); + problem.options.simdLen=simdLen; + + typename CRHMCWalk::template parameters + < + NT, + NegativeGradientFunctor + > params(F, P.dimension(), problem.options); + + + if (input.df.params.eta > 0) { + params.eta = input.df.params.eta; + } + + + 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, @@ -450,28 +436,32 @@ void exponential_sampling(PointList &randPoints, const unsigned int &rnum, const Point &c, const NT &a, - const Point &starting_point, + Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk < - Polytope, + Polytope, RandomNumberGenerator > 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); } @@ -494,28 +484,25 @@ void exponential_sampling(PointList &randPoints, const Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator > walk; + + walk walkInstance(P, starting_point, rng, WalkType.param); - 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, @@ -533,28 +520,26 @@ void exponential_sampling(PointList &randPoints, const Point &c, const NT &a, const NT &eta, - const Point &starting_point, + Point &starting_point, unsigned int const& nburns) { - 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); } @@ -578,25 +563,23 @@ void exponential_sampling(PointList &randPoints, const Point &starting_point, unsigned int const& nburns) { - - typedef typename WalkTypePolicy::template Walk + + typename WalkTypePolicy::template Walk < Polytope, RandomNumberGenerator - > walk; - - PushBackWalkPolicy push_back_policy; + > 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); }