From cfa7cf8df5eeef894b721482a3ed6b30eba667ca Mon Sep 17 00:00:00 2001 From: vgnecula Date: Fri, 16 Aug 2024 14:37:27 -0400 Subject: [PATCH] Improved CRHMC burnIn logic, with additional random hpoly examples --- .../volume_example.cpp | 12 ++ ...e_cooling_nonspherical_gaussians_crhmc.hpp | 152 ++++++++++++------ 2 files changed, 113 insertions(+), 51 deletions(-) diff --git a/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp b/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp index 8789abee8..c7906b5d4 100644 --- a/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp +++ b/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp @@ -26,6 +26,7 @@ typedef double NT; typedef Cartesian Kernel; typedef typename Kernel::Point Point; typedef BoostRandomNumberGenerator RandomNumberGenerator; +typedef boost::mt19937 PolyRNGType; typedef HPolytope HPOLYTOPE; NT nonspherical_crhmc_volume(HPOLYTOPE& polytope) { @@ -65,5 +66,16 @@ int main() { std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(skinnycube3) << "\n"; std::cout << "Expected Volume: " << 200 * std::pow(2, 2) << "\n\n"; + HPOLYTOPE P3 = random_hpoly(3, 12, false); + std::cout << "Random 3D Hpoly \n"; + std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(P3) << "\n"; + std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(P3) << "\n"; + std::cout << "Expected Volume: " << "N/A" << "\n\n"; + + HPOLYTOPE P4 = random_hpoly(4, 16, false); + std::cout << "Random 4D Hpoly \n"; + std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(P4) << "\n"; + std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(P4) << "\n"; + std::cout << "Expected Volume: " << "N/A" << "\n\n"; return 0; } \ No newline at end of file diff --git a/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp b/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp index 47978ed06..21e65cef0 100644 --- a/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp +++ b/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp @@ -24,6 +24,7 @@ // Gaussian Anealling // Compute a_{i+1} when a_i is given + template < typename CRHMCWalkType, @@ -32,6 +33,36 @@ template typename Grad, typename Func, typename CrhmcProblem, + typename Point, + typename NT, + typename RandomNumberGenerator +> +void burnIn(Point &p, + const unsigned int& walk_length, + RandomNumberGenerator& rng, + Grad& g, + Func& f, + crhmc_walk_params& parameters, + CrhmcProblem& problem, + CRHMCWalkType& crhmc_walk, + NT burnIn_sample) +{ + std::list randPoints; + + // burnIn + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + CRHMCRandomPointGenerator::apply(problem, p, burnIn_sample, walk_length, randPoints, + push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); + +} + + +template +< + int simdLen, typename Polytope, typename Point, typename NT, @@ -39,20 +70,74 @@ template typename RandomNumberGenerator > NT get_next_gaussian(Polytope& P, - Point &p, + Point& p, NT const& a, const unsigned int &N, const NT &ratio, const NT &C, const unsigned int& walk_length, RandomNumberGenerator& rng, - Grad& g, - Func& f, - crhmc_walk_params& parameters, - CrhmcProblem& problem, - CRHMCWalkType& crhmc_walk, MT const& inv_covariance_matrix) { + typedef typename NonSphericalGaussianFunctor::FunctionFunctor Func; + typedef typename NonSphericalGaussianFunctor::GradientFunctor Grad; + typedef typename NonSphericalGaussianFunctor::HessianFunctor Hess; + typedef typename NonSphericalGaussianFunctor::parameters func_params; + + typedef crhmc_input> Input; + //typedef crhmc_input Input; + + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename CRHMCWalk::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename CRHMCWalk::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + // Create the CRHMC problem for this variance + int dimension = P.dimension(); + func_params f_params = func_params(Point(dimension), a, -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + Hess h(f_params); + ZeroFunctor zerof; + + Input input = convert2crhmc_input>(P, f, g, zerof); + + //Input input = convert2crhmc_input(P, f, g, h); + + typedef crhmc_problem CrhmcProblem; + + CrhmcProblem problem = CrhmcProblem(input); + + if(problem.terminate) { return 0;} + + problem.options.simdLen = simdLen; + + crhmc_walk_params params(input.df, p.dimension(), problem.options); + + // Create the walk object for this problem + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + burnIn( + p, walk_length, rng, g, f, params, problem, walk, 2.5 * N); + NT last_a = a; NT last_ratio = 0.1; NT k = 1.0; @@ -67,7 +152,7 @@ NT get_next_gaussian(Polytope& P, typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints, - push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); + push_back_policy, rng, g, f, params, walk, simdLen, raw_output); while (!done) { NT new_a = last_a * std::pow(ratio, k); @@ -96,41 +181,6 @@ NT get_next_gaussian(Polytope& P, return last_a * std::pow(ratio, k); } -template -< - typename CRHMCWalkType, - typename crhmc_walk_params, - int simdLen, - typename Grad, - typename Func, - typename CrhmcProblem, - typename Point, - typename NT, - typename RandomNumberGenerator -> -void burnIn(Point &p, - const unsigned int& walk_length, - RandomNumberGenerator& rng, - Grad& g, - Func& f, - crhmc_walk_params& parameters, - CrhmcProblem& problem, - CRHMCWalkType& crhmc_walk, - NT burnIn_sample) -{ - std::list randPoints; - - // burnIn - PushBackWalkPolicy push_back_policy; - bool raw_output = false; - typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; - - CRHMCRandomPointGenerator::apply(problem, p, burnIn_sample, walk_length, randPoints, - push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); - -} - - // Compute the sequence of non spherical gaussians template< int simdLen, @@ -188,7 +238,6 @@ void compute_annealing_schedule(Polytope Pin_copy, auto ball = P_copy.ComputeInnerBall(); P_copy.shift(ball.first.getCoefficients()); // when using max_ellipsoid for rounding this center is the origin, but if we use other covariances this is different than the origin get_first_gaussian(P_copy, frac, ball.second, error, a_vals); // the function included from volume_cooling_gaussians.hpp (spherical gaussians) - #ifdef VOLESTI_DEBUG std::cout << "first gaussian computed " << a_vals[0] << std::endl; #endif @@ -213,9 +262,9 @@ void compute_annealing_schedule(Polytope Pin_copy, Func initial_f(initial_f_params); Grad initial_g(initial_f_params); Hess initial_h(initial_f_params); - ZeroFunctor zerof; + ZeroFunctor initial_zerof; - Input initial_input = convert2crhmc_input>(Pin_copy, initial_f, initial_g, zerof); + Input initial_input = convert2crhmc_input>(Pin_copy, initial_f, initial_g, initial_zerof); CrhmcProblem initial_problem = CrhmcProblem(initial_input); Point initial_p = Point(initial_problem.center); @@ -234,7 +283,11 @@ void compute_annealing_schedule(Polytope Pin_copy, while (true) { - + // Compute the next gaussian + NT next_a = get_next_gaussian( + Pin_copy, start_point, a_vals[it], N, ratio, C, walk_length, rng, inv_covariance_matrix); + + // Decide if this is the last one NT curr_fn = 0; NT curr_its = 0; auto steps = totalSteps; @@ -248,6 +301,7 @@ void compute_annealing_schedule(Polytope Pin_copy, Func f(f_params); Grad g(f_params); Hess h(f_params); + ZeroFunctor zerof; Input input = convert2crhmc_input>(Pin_copy, f, g, zerof); @@ -270,11 +324,6 @@ void compute_annealing_schedule(Polytope Pin_copy, // Create the walk object for this problem CRHMCWalkType walk = CRHMCWalkType(problem, start_point, input.df, input.f, params); - // Compute the next gaussian - NT next_a = get_next_gaussian( - Pin_copy, start_point, a_vals[it], N, ratio, C, walk_length, rng, g, f, params, problem, walk, inv_covariance_matrix); - - // Compute some ratios to decide if this is the last gaussian for (unsigned int j = 0; j < steps; j++) { @@ -400,6 +449,7 @@ double non_spherical_crhmc_volume_cooling_gaussians(Polytope& Pin, auto L = lltOfA.matrixL(); P.linear_transformIt(L); + // Initialize the gaussian_annealing_parameters struct non_gaussian_annealing_parameters parameters(P.dimension());