From 2a2ec8e36a9d93a8002435f687e1fd507dfc4c70 Mon Sep 17 00:00:00 2001 From: vgnecula Date: Sun, 30 Jun 2024 22:38:04 -0400 Subject: [PATCH] CRHMC Volume Cooling Gaussians --- examples/crhmc_sampling/CMakeLists.txt | 2 + examples/crhmc_sampling/volume_example.cpp | 40 ++ .../volume/volume_cooling_gaussians_crhmc.hpp | 616 ++++++++++++++++++ 3 files changed, 658 insertions(+) create mode 100644 examples/crhmc_sampling/volume_example.cpp create mode 100644 include/volume/volume_cooling_gaussians_crhmc.hpp diff --git a/examples/crhmc_sampling/CMakeLists.txt b/examples/crhmc_sampling/CMakeLists.txt index 91fe46a7f..6dc810179 100644 --- a/examples/crhmc_sampling/CMakeLists.txt +++ b/examples/crhmc_sampling/CMakeLists.txt @@ -125,5 +125,7 @@ else () TARGET_LINK_LIBRARIES(sampling_functions QD_LIB ${MKL_LINK} ${LP_SOLVE}) add_executable (simple_crhmc simple_crhmc.cpp) TARGET_LINK_LIBRARIES(simple_crhmc QD_LIB ${MKL_LINK} ${LP_SOLVE}) + add_executable(volume_example volume_example.cpp) + target_link_libraries(volume_example QD_LIB ${MKL_LINK} ${LP_SOLVE}) endif() diff --git a/examples/crhmc_sampling/volume_example.cpp b/examples/crhmc_sampling/volume_example.cpp new file mode 100644 index 000000000..22e1ccd56 --- /dev/null +++ b/examples/crhmc_sampling/volume_example.cpp @@ -0,0 +1,40 @@ +#include "Eigen/Eigen" +#include +#include "cartesian_geom/cartesian_kernel.h" +#include "convex_bodies/hpolytope.h" +#include "generators/known_polytope_generators.h" + +#include "random_walks/random_walks.hpp" +#include "volume/volume_sequence_of_balls.hpp" +#include "volume/volume_cooling_gaussians_crhmc.hpp" +#include "volume/volume_cooling_balls.hpp" + +#include +#include +#include "misc/misc.h" + +typedef double NT; +typedef Cartesian Kernel; +typedef typename Kernel::Point Point; +typedef BoostRandomNumberGenerator RNGType; +typedef HPolytope HPOLYTOPE; + +template +void calculateAndVerifyVolume(HPOLYTOPE& polytope) { + int walk_len = 100; + NT e = 0.1; + + NT volume = volume_cooling_gaussians(polytope, e, walk_len); + + std::cout << "Volume " << volume << std::endl; +} + +int main() { + + HPOLYTOPE HP = generate_simplex(2,false); + std::cout << "HPoly: " << std::endl; + HP.print(); + calculateAndVerifyVolume<4>(HP); + + return 0; +} \ No newline at end of file diff --git a/include/volume/volume_cooling_gaussians_crhmc.hpp b/include/volume/volume_cooling_gaussians_crhmc.hpp new file mode 100644 index 000000000..6a48f49fa --- /dev/null +++ b/include/volume/volume_cooling_gaussians_crhmc.hpp @@ -0,0 +1,616 @@ +#ifndef VOLUME_COOLING_GAUSSIANS_HPP +#define VOLUME_COOLING_GAUSSIANS_HPP + +//#define VOLESTI_DEBUG + +#include +#include +#include +#include +#include + +#include "cartesian_geom/cartesian_kernel.h" +#include "random_walks/gaussian_helpers.hpp" +#include "random_walks/gaussian_ball_walk.hpp" +#include "random_walks/gaussian_cdhr_walk.hpp" +#include "sampling/random_point_generators.hpp" +#include "volume/math_helpers.hpp" + +//new include for crhmc +#include "Eigen/Eigen" +#include "cartesian_geom/cartesian_kernel.h" +#include "volume/sampling_policies.hpp" +#include "ode_solvers/ode_solvers.hpp" +#include "preprocess/crhmc/crhmc_input.h" +#include "preprocess/crhmc/crhmc_problem.h" +#include "sampling/random_point_generators.hpp" +#include "sampling/sampling.hpp" +#include "random.hpp" +#include +#include "random_walks/random_walks.hpp" +#include "generators/known_polytope_generators.h" + + +// define types +using NT = double; +using Kernel = Cartesian; +using Point = typename Kernel::Point; +using Func = GaussianFunctor::FunctionFunctor; +using Grad = GaussianFunctor::GradientFunctor; +using Hess = GaussianFunctor::HessianFunctor; + +using NegativeLogprobFunctor = GaussianFunctor::FunctionFunctor; //Func +using NegativeGradientFunctor = GaussianFunctor::GradientFunctor; //Grad +using HessianFunctor = GaussianFunctor::HessianFunctor; //Hess + +using PolytopeType = HPolytope; +using MT = PolytopeType::MT; +using VT = Eigen::Matrix; +using func_params = GaussianFunctor::parameters; +using RNG = BoostRandomNumberGenerator; + +//Param building -> see sampling_functions.cpp +using Input = crhmc_input; +using CrhmcProblem = crhmc_problem; + + +/////////////////// Helpers for random walks + +template +struct update_delta +{ + template + static void apply(WalkType& walk, NT delta) {} +}; + +template +struct update_delta> +{ + template + static void apply(GaussianBallWalk::Walk walk, + NT delta) + { + walk.update_delta(delta); + } +}; + + +////////////////////////////// Algorithms + +// Gaussian Anealling + +// Implementation is based on algorithm from paper "A practical volume algorithm", +// Springer-Verlag Berlin Heidelberg and The Mathematical Programming Society 2015 +// Ben Cousins, Santosh Vempala + +// Compute the first variance a_0 for the starting gaussian +template +void get_first_gaussian(Polytope& P, + NT const& frac, + NT const& chebychev_radius, + NT const& error, + std::vector & a_vals) +{ + // if tol is smaller than 1e-6 no convergence can be obtained when float is used + NT tol = std::is_same::value ? 0.001 : 0.0000001; + + std::vector dists = P.get_dists(chebychev_radius); + NT lower = 0.0; + NT upper = 1.0; + + // Compute an upper bound for a_0 + unsigned int i; + const unsigned int maxiter = 10000; + for (i= 1; i <= maxiter; ++i) { + NT sum = 0.0; + for (auto it = dists.begin(); it != dists.end(); ++it) + { + sum += std::exp(-upper * std::pow(*it, 2.0)) + / (2.0 * (*it) * std::sqrt(M_PI * upper)); + } + + if (sum > frac * error) + { + upper = upper * 10; + } else { + break; + } + } + + if (i == maxiter) { +#ifdef VOLESTI_DEBUG + std::cout << "Cannot obtain sharp enough starting Gaussian" << std::endl; +#endif + return; + } + + //get a_0 with binary search + while (upper - lower > tol) + { + NT mid = (upper + lower) / 2.0; + NT sum = 0.0; + for (auto it = dists.begin(); it != dists.end(); ++it) { + sum += std::exp(-mid * std::pow(*it, 2.0)) + / (2.0 * (*it) * std::sqrt(M_PI * mid)); + } + + if (sum < frac * error) { + upper = mid; + } else { + lower = mid; + } + } + a_vals.push_back((upper + lower) / NT(2.0)); +} + + +// Compute a_{i+1} when a_i is given +template +< + typename WalkType, + typename walk_params, + typename RandomPointGenerator, + int simdLen, + typename Polytope, + typename Point, + typename NT, + typename RandomNumberGenerator +> +NT get_next_gaussian(Polytope& 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, + walk_params& parameters, + CrhmcProblem& problem, + WalkType& crhmc_walk) +{ + + NT last_a = a; + NT last_ratio = 0.1; + //k is needed for the computation of the next variance a_{i+1} = a_i * (1-1/d)^k + NT k = 1.0; + const NT tol = 0.00001; + bool done=false; + std::vector fn(N,NT(0.0)); + std::list randPoints; + typedef typename std::vector::iterator viterator; + + //sample N points + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + + RandomPointGenerator::apply(problem, p, N, walk_length, randPoints, + push_back_policy, rng, g, f, parameters, crhmc_walk, simdLen, raw_output); + + while (!done) + { + NT new_a = last_a * std::pow(ratio,k); + + auto fnit = fn.begin(); + for (auto pit=randPoints.begin(); pit!=randPoints.end(); ++pit, fnit++) + { + *fnit = eval_exp(*pit, new_a)/eval_exp(*pit, last_a); + } + std::pair mv = get_mean_variance(fn); + + // Compute a_{i+1} + if (mv.second/(mv.first * mv.first)>=C || mv.first/last_ratio<1.0+tol) + { + if (k != 1.0) + { + k = k / 2; + } + done = true; + } else { + k = 2 * k; + } + last_ratio = mv.first; + } + return last_a * std::pow(ratio, k); +} + +// Compute the sequence of spherical gaussians +template +< + typename WalkType, + typename walk_params, + typename RandomPointGenerator, + int simdLen, + typename Polytope, + typename NT, + typename RandomNumberGenerator +> +void compute_annealing_schedule(Polytope& P, + NT const& ratio, + NT const& C, + NT const& frac, + unsigned int const& N, + unsigned int const& walk_length, + NT const& chebychev_radius, + NT const& error, + std::vector& a_vals, + RandomNumberGenerator& rng) +{ + + + typedef typename Polytope::PointType Point; + typedef typename Polytope::VT VT; + + // Compute the first gaussian + get_first_gaussian(P, frac, chebychev_radius, error, a_vals); + NT a_stop = 0.0; + const NT tol = 0.001; + unsigned int it = 0; + unsigned int n = P.dimension(); + const unsigned int totalSteps = ((int)150/((1.0 - frac) * error))+1; + + if (a_vals[0](P, f, g, h); + + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + + Point p = Point(problem.center); + + if(problem.terminate){return;} + + 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; + } + + int dim; + dim = p.dimension(); + + //create the walk object for this problem + WalkType walk = WalkType(problem, p, input.df, input.f, params); + + //TODO: test update delta here? + update_delta + ::apply(walk, 4.0 * chebychev_radius + / std::sqrt(std::max(NT(1.0), a_vals[it]) * NT(n))); + + // Compute the next gaussian + NT next_a = get_next_gaussian + (P, p, a_vals[it], N, ratio, C, walk_length, rng, g, f, params, problem, walk); + +#ifdef VOLESTI_DEBUG + std::cout<<"Next Gaussian " << next_a <0 && curr_fn/curr_its>(1.0+tol)) + { + a_vals.push_back(next_a); + it++; + } else if (next_a <= 0) + { + a_vals.push_back(a_stop); + it++; + break; + } else { + a_vals[it] = a_stop; + break; + } + } +#ifdef VOLESTI_DEBUG + std::cout<<"first gaussian after WHILE "<< a_vals[0] < +struct gaussian_annealing_parameters +{ + gaussian_annealing_parameters(unsigned int d) + : frac(0.1) + , ratio(NT(1)-NT(1)/(NT(d))) + , C(NT(2)) + , N(500 * ((int) C) + ((int) (d * d / 2))) + , W(6*d*d+800) + {} + + NT frac; + NT ratio; + NT C; + unsigned int N; + unsigned int W; +}; + +template +< + typename WalkTypePolicy, + typename Polytope, + typename RandomNumberGenerator, + int simdLen +> +double volume_cooling_gaussians(Polytope& Pin, + RandomNumberGenerator& rng, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + using Solver = + ImplicitMidpointODESolver; + + typedef typename WalkTypePolicy::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + NegativeGradientFunctor, + NegativeLogprobFunctor, + Solver + > WalkType; + + typedef typename WalkTypePolicy::template parameters + < + NT, + NegativeGradientFunctor + > walk_params; + + typedef CrhmcRandomPointGenerator RandomPointGenerator; + + //const NT maxNT = std::numeric_limits::max();//1.79769e+308; + //const NT minNT = std::numeric_limits::min();//-1.79769e+308; + + auto P(Pin); //copy and work with P because we are going to shift + unsigned int n = P.dimension(); + unsigned int m = P.num_of_hyperplanes(); + gaussian_annealing_parameters parameters(P.dimension()); + //RandomNumberGenerator rng(n); + + // Consider Chebychev center as an internal point + auto InnerBall = P.ComputeInnerBall(); + if (InnerBall.second < 0.0) return -1.0; + + Point c = InnerBall.first; + NT radius = InnerBall.second; + + // Move the chebychev center to the origin and apply the same shifting to the polytope + P.shift(c.getCoefficients()); + + // Computing the sequence of gaussians +#ifdef VOLESTI_DEBUG + std::cout<<"\n\nComputing annealing...\n"< a_vals; + NT ratio = parameters.ratio; + NT C = parameters.C; + unsigned int N = parameters.N; + + compute_annealing_schedule + < + WalkType, + walk_params, + RandomPointGenerator, + simdLen + >(P, ratio, C, parameters.frac, N, walk_length, radius, error, a_vals, rng); + +#ifdef VOLESTI_DEBUG + std::cout<<"All the variances of schedule_annealing computed in = " + << (double)clock()/(double)CLOCKS_PER_SEC-tstart2<<" sec"< last_W2(W,0); + std::vector fn(mm,0); + std::vector its(mm,0); + VT lamdas; + lamdas.setZero(m); + NT vol = std::pow(M_PI/a_vals[0], (NT(n))/2.0); + unsigned int i=0; + + typedef typename std::vector::iterator viterator; + viterator itsIt = its.begin(); + viterator avalsIt = a_vals.begin(); + viterator minmaxIt; + + +#ifdef VOLESTI_DEBUG + std::cout<<"volume of the first gaussian = "<::min(); + NT max_val = std::numeric_limits::max(); + unsigned int min_index = W-1; + unsigned int max_index = W-1; + unsigned int index = 0; + unsigned int min_steps = 0; + std::vector last_W = last_W2; + + // Set the radius for the ball walk + //creating the walk object + int dimension = P.dimension(); + func_params f_params = func_params(Point(dimension), *avalsIt, 1); + + Func f(f_params); + Grad g(f_params); + Hess h(f_params); + + //create the crhmc problem + Input input = convert2crhmc_input(P, f, g, h); + + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + + Point p = Point(problem.center); + + if(problem.terminate){return 0;} + + 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; + } + + WalkType walk = WalkType(problem, p, input.df, input.f, params); + + + update_delta + ::apply(walk, 4.0 * radius + / std::sqrt(std::max(NT(1.0), *avalsIt) * NT(n))); + + while (!done || (*itsIt)= max_val) + { + max_val = val; + max_index = index; + } else if (max_index == index) + { + minmaxIt = std::max_element(last_W.begin(), last_W.end()); + max_val = *minmaxIt; + max_index = std::distance(last_W.begin(), minmaxIt); + } + + if ( (max_val-min_val)/max_val <= curr_eps/2.0 ) + { + done=true; + } + + index = index%W + 1; + if (index == W) index = 0; + } +#ifdef VOLESTI_DEBUG + std::cout << "ratio " << i << " = " << (*fnIt) / (*itsIt) + << " N_" << i << " = " << *itsIt << std::endl; +#endif + vol *= ((*fnIt) / (*itsIt)); + } + +#ifdef VOLESTI_DEBUG + NT sum_of_steps = 0.0; + for(viterator it = its.begin(); it != its.end(); ++it) { + sum_of_steps += *it; + } + auto steps= int(sum_of_steps); + std::cout<<"\nTotal number of steps = "<, + typename Polytope, + int simdLen +> +double volume_cooling_gaussians(Polytope &Pin, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + RandomNumberGenerator rng(Pin.dimension()); + return volume_cooling_gaussians(Pin, rng, error, walk_length); +} + + +template +< + typename WalkTypePolicy = GaussianCDHRWalk, + typename RandomNumberGenerator = BoostRandomNumberGenerator, + typename Polytope, + int simdLen +> +double volume_cooling_gaussians(Polytope &Pin, + Cartesian::Point const& interior_point, + unsigned int const& walk_length = 1, + double const& error = 0.1) +{ + RandomNumberGenerator rng(Pin.dimension()); + Pin.set_interior_point(interior_point); + + return volume_cooling_gaussians(Pin, rng, error, walk_length); +} + +#endif // VOLUME_COOLING_GAUSSIANS_HPP \ No newline at end of file