From e19ce969c94add52965eac5698c41d939ad6c86c Mon Sep 17 00:00:00 2001 From: vgnecula Date: Sat, 20 Apr 2024 01:24:58 -0400 Subject: [PATCH] BurnIn Handling, same style --- include/sampling/sampling.hpp | 368 +++++++++++++++------------------- 1 file changed, 163 insertions(+), 205 deletions(-) diff --git a/include/sampling/sampling.hpp b/include/sampling/sampling.hpp index 696186098..5b16e892c 100644 --- a/include/sampling/sampling.hpp +++ b/include/sampling/sampling.hpp @@ -33,31 +33,24 @@ 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 - < - 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 < @@ -73,27 +66,24 @@ 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 - < - 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); } @@ -113,27 +103,22 @@ void uniform_sampling_boundary(PointList &randPoints, const Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk - < - Polytope, - RandomNumberGenerator - > walk; - - //RandomNumberGenerator rng(P.dimension()); - PushBackWalkPolicy push_back_policy; + typedef typename WalkTypePolicy::template Walk Walk; - 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/2; ++i) { + walk.apply(P, current_point, dummy_point, walk_len, rng); + randPoints.push_back(current_point); + current_point = dummy_point; + } } @@ -156,27 +141,20 @@ void gaussian_sampling(PointList &randPoints, 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); + } } @@ -199,25 +177,20 @@ void gaussian_sampling(PointList &randPoints, 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 < @@ -241,50 +214,33 @@ void logconcave_sampling(PointList &randPoints, NegativeGradientFunctor &F, NegativeLogprobFunctor &f) { - typedef typename WalkTypePolicy::template Walk - < - Point, - Polytope, - RandomNumberGenerator, - NegativeGradientFunctor, - NegativeLogprobFunctor, - Solver - > walk; - - typedef typename WalkTypePolicy::template parameters - < - NT, - NegativeGradientFunctor - > walk_params; + + typename WalkTypePolicy::template Walk< + Polytope, RandomNumberGenerator, NegativeGradientFunctor, NegativeLogprobFunctor, Solver + > 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 +265,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 +411,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 +459,25 @@ void exponential_sampling(PointList &randPoints, const Point &starting_point, unsigned int const& nburns) { - typedef typename WalkTypePolicy::template Walk < - Polytope, + 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 +495,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, + 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 +538,23 @@ void exponential_sampling(PointList &randPoints, const Point &starting_point, unsigned int const& nburns) { - - typedef typename WalkTypePolicy::template Walk + + typename WalkTypePolicy::template Walk < - Polytope, + 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); }