diff --git a/examples/crhmc_cooling_nonspherical_gaussians/CMakeLists.txt b/examples/crhmc_cooling_nonspherical_gaussians/CMakeLists.txt new file mode 100644 index 000000000..d3577f9c2 --- /dev/null +++ b/examples/crhmc_cooling_nonspherical_gaussians/CMakeLists.txt @@ -0,0 +1,132 @@ +# VolEsti (volume computation and sampling library) + +# Copyright (c) 2012-2024 Vissarion Fisikopoulos +# Copyright (c) 2018-2024 Apostolos Chalkis +# Copyright (c) 2024 Vladimir Necula + +# Contributed and/or modified by Vladimir Necula, as part of Google Summer of +# Code 2024 program. + +# Licensed under GNU LGPL.3, see LICENCE file + +project( VolEsti ) + + +CMAKE_MINIMUM_REQUIRED(VERSION 3.11) + +set(CMAKE_ALLOW_LOOSE_LOOP_CONSTRUCTS true) + +# Locate Intel MKL root (in case it is enabled) +if (APPLE) + set(MKLROOT /opt/intel/oneapi/mkl/latest) +elseif(UNIX) + #set(MKLROOT /opt/intel/oneapi/mkl/latest) + set(MKLROOT $ENV{HOME}/intel/mkl) +endif() + +if(COMMAND cmake_policy) + cmake_policy(SET CMP0003 NEW) +endif(COMMAND cmake_policy) + + +option(DISABLE_NLP_ORACLES "Disable non-linear oracles (used in collocation)" ON) +option(BUILTIN_EIGEN "Use eigen from ../../external" OFF) +option(USE_MKL "Use MKL library to build eigen" OFF) + + +if(DISABLE_NLP_ORACLES) + add_definitions(-DDISABLE_NLP_ORACLES) +else() + find_library(IFOPT NAMES libifopt_core.so PATHS /usr/local/lib) + find_library(IFOPT_IPOPT NAMES libifopt_ipopt.so PATHS /usr/local/lib) + find_library(GMP NAMES libgmp.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu) + find_library(MPSOLVE NAMES libmps.so PATHS /usr/local/lib) + find_library(PTHREAD NAMES libpthread.so PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu) + find_library(FFTW3 NAMES libfftw3.so.3 PATHS /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu) + + if (NOT IFOPT) + + message(FATAL_ERROR "This program requires the ifopt library, and will not be compiled.") + + elseif (NOT GMP) + + message(FATAL_ERROR "This program requires the gmp library, and will not be compiled.") + + elseif (NOT MPSOLVE) + + message(FATAL_ERROR "This program requires the mpsolve library, and will not be compiled.") + + elseif (NOT FFTW3) + + message(FATAL_ERROR "This program requires the fftw3 library, and will not be compiled.") + + else() + message(STATUS "Library ifopt found: ${IFOPT}") + message(STATUS "Library gmp found: ${GMP}") + message(STATUS "Library mpsolve found: ${MPSOLVE}") + message(STATUS "Library fftw3 found:" ${FFTW3}) + + endif(NOT IFOPT) + +endif(DISABLE_NLP_ORACLES) + +include("../../external/cmake-files/Eigen.cmake") +GetEigen() + +include("../../external/cmake-files/Boost.cmake") +GetBoost() + +include("../../external/cmake-files/LPSolve.cmake") +GetLPSolve() + +include("../../external/cmake-files/QD.cmake") +GetQD() + +# Find lpsolve library +find_library(LP_SOLVE NAMES liblpsolve55.so PATHS /usr/lib/lp_solve) + +if (NOT LP_SOLVE) + message(FATAL_ERROR "This program requires the lp_solve library, and will not be compiled.") +else () + message(STATUS "Library lp_solve found: ${LP_SOLVE}") + + set(CMAKE_EXPORT_COMPILE_COMMANDS "ON") + + if (USE_MKL) + find_library(BLAS NAMES libblas.so libblas.dylib PATHS /usr/local/Cellar/lapack/3.9.1_1/lib /usr/lib/x86_64-linux-gnu /usr/lib/i386-linux-gnu /usr/local/Cellar/openblas/0.3.15_1/lib /usr/lib) + find_library(GFORTRAN NAME libgfortran.dylib PATHS /usr/local/Cellar/gcc/10.2.0_4/lib/gcc/10) + find_library(LAPACK NAME liblapack.dylib PATHS /usr/lib) + find_library(OPENMP NAME libiomp5.dylib PATHS /opt/intel/oneapi/compiler/2021.1.1/mac/compiler/lib) + + include_directories (BEFORE ${MKLROOT}/include) + set(PROJECT_LIBS ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} ${GFORTRAN_LIBRARIES}) + set(MKL_LINK "-L${MKLROOT}/lib -Wl,-rpath,${MKLROOT}/lib -lmkl_intel_ilp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl") + add_definitions(-DEIGEN_USE_MKL_ALL) + else() + set(MKL_LINK "") + endif(USE_MKL) + + include_directories (BEFORE ../../external) + include_directories (BEFORE ../../external/minimum_ellipsoid) + include_directories (BEFORE ../../include/) + + # for Eigen + if (${CMAKE_VERSION} VERSION_LESS "3.12.0") + add_compile_options(-D "EIGEN_NO_DEBUG") + else () + add_compile_definitions("EIGEN_NO_DEBUG") + endif () + + + add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17") # enable C++17 standard + set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING") + add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING" ${ADDITIONAL_FLAGS}) # optimization of the compiler + #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl") + add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm") + add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl") + add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-DBOOST_NO_AUTO_PTR") + + add_executable(volume_example volume_example.cpp) + target_link_libraries(volume_example QD_LIB ${MKL_LINK} ${LP_SOLVE}) + +endif() \ No newline at end of file diff --git a/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp b/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp new file mode 100644 index 000000000..c7906b5d4 --- /dev/null +++ b/examples/crhmc_cooling_nonspherical_gaussians/volume_example.cpp @@ -0,0 +1,81 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2024 Vissarion Fisikopoulos +// Copyright (c) 2018-2024 Apostolos Chalkis +// Copyright (c) 2024 Vladimir Necula + +// Contributed and/or modified by Vladimir Necula, as part of Google Summer of +// Code 2024 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +#include "generators/known_polytope_generators.h" +#include "generators/h_polytopes_generator.h" +#include "random_walks/random_walks.hpp" +#include "volume/volume_cooling_nonspherical_gaussians_crhmc.hpp" +#include "volume/volume_cooling_gaussians.hpp" +#include +#include +#include +#include +#include "misc/misc.h" + +const unsigned int FIXED_SEED = 42; + +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) { + int walk_len = 10; + NT e = 0.1; + RandomNumberGenerator rng; + // NT volume = volume_cooling_gaussians(polytope, e, walk_len); + NT volume = non_spherical_crhmc_volume_cooling_gaussians(polytope, rng, e, walk_len); + return volume; +} + +NT spherical_gaussians_volume(HPOLYTOPE& polytope) { + int walk_len = 10; + NT e = 0.1; + RandomNumberGenerator rng; + NT volume = volume_cooling_gaussians(polytope, e, walk_len); + return volume; +} + +int main() { + + HPOLYTOPE cube3 = generate_cube(3, false); + std::cout << "Cube3 \n"; + std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(cube3) << "\n"; + std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(cube3) << "\n"; + std::cout << "Expected Volume: " << std::pow(2, 3) << "\n\n"; + + HPOLYTOPE cube4 = generate_cube(4, false); + std::cout << "Cube4 \n"; + std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(cube4) << "\n"; + std::cout << "Calculated Volume With CRHMC: " << nonspherical_crhmc_volume(cube4) << "\n"; + std::cout << "Expected Volume: " << std::pow(2, 4) << "\n\n"; + + HPOLYTOPE skinnycube3 = generate_skinny_cube(3, false); + std::cout << "SkinnyCube3 \n"; + std::cout << "Calculated Volume With Gaussian CDHR: " << spherical_gaussians_volume(skinnycube3) << "\n"; + 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/ode_solvers/implicit_midpoint.hpp b/include/ode_solvers/implicit_midpoint.hpp index 786734d84..7f13bfe9b 100644 --- a/include/ode_solvers/implicit_midpoint.hpp +++ b/include/ode_solvers/implicit_midpoint.hpp @@ -172,6 +172,8 @@ struct ImplicitMidpointODESolver { stream << '\n'; } } + + NT get_eta() const { return eta; } }; #endif diff --git a/include/ode_solvers/oracle_functors.hpp b/include/ode_solvers/oracle_functors.hpp index 35fc8887b..ab4546b69 100644 --- a/include/ode_solvers/oracle_functors.hpp +++ b/include/ode_solvers/oracle_functors.hpp @@ -394,4 +394,77 @@ struct HessianFunctor { }; +struct NonSphericalGaussianFunctor { + template + struct parameters { + Point x0; + NT a; + NT eta; + unsigned int order; + NT L; + NT m; + NT kappa; + Eigen::Matrix inv_covariance_matrix; + + parameters(Point x0_, NT a_, NT eta_, Eigen::Matrix inv_covariance_matrix_) + : x0(x0_), a(a_), eta(eta_), order(2), + L(a_), m(a_), kappa(1), + inv_covariance_matrix(inv_covariance_matrix_) {}; + }; + + template + struct GradientFunctor { + typedef typename Point::FT NT; + typedef std::vector pts; + + parameters& params; + + GradientFunctor(parameters& params_) : params(params_) {}; + + Point operator()(unsigned int const& i, pts const& xs, NT const& t) const { + if (i == params.order - 1) { + Point y = xs[0] - params.x0; + Eigen::Matrix result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients(); + return Point(result); + } else { + return xs[i + 1]; + } + } + Point operator()(Point const& x) { + Point y = x - params.x0; + Eigen::Matrix result = -2.0 * params.a * params.inv_covariance_matrix * y.getCoefficients(); + return Point(result); + } + }; + + template + struct FunctionFunctor { + typedef typename Point::FT NT; + + parameters& params; + + FunctionFunctor(parameters& params_) : params(params_) {}; + + NT operator()(Point const& x) const { + Point y = x - params.x0; + Eigen::Matrix y_coeffs = y.getCoefficients(); + return params.a * y_coeffs.dot(params.inv_covariance_matrix * y_coeffs); + } + }; + + template + struct HessianFunctor { + typedef typename Point::FT NT; + + parameters& params; + + HessianFunctor(parameters& params_) : params(params_) {}; + + Point operator()(Point const& x) const { + Eigen::Matrix result = 2.0 * params.a * params.inv_covariance_matrix; + return Point(result); + } + }; +}; + #endif diff --git a/include/preprocess/crhmc/opts.h b/include/preprocess/crhmc/opts.h index 573a69dbd..fffd29c58 100644 --- a/include/preprocess/crhmc/opts.h +++ b/include/preprocess/crhmc/opts.h @@ -44,6 +44,7 @@ template class opts { bool DynamicWeight = true; //Enable the use of dynamic weights for each variable when sampling bool DynamicStepSize = true; // Enable adaptive step size that avoids low acceptance probability bool DynamicRegularizer = true; //Enable the addition of a regularization term + bool etaInitialize = true; // enable eta initialization Type regularization_factor=1e-20; /*Dynamic step choices*/ Type warmUpStep = 10; diff --git a/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp b/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp index 367211962..3daee1c34 100644 --- a/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp +++ b/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp @@ -50,7 +50,7 @@ class dynamic_step_size { consecutiveBadStep = IVT::Zero(simdLen); rejectSinceShrink = VT::Zero(simdLen); - if (options.warmUpStep > 0) { + if (options.warmUpStep > 0 && options.etaInitialize) { eta = 1e-3; } else { warmupFinished = true; diff --git a/include/random_walks/crhmc/crhmc_walk.hpp b/include/random_walks/crhmc/crhmc_walk.hpp index 21ed8d6bb..9f1dd24af 100644 --- a/include/random_walks/crhmc/crhmc_walk.hpp +++ b/include/random_walks/crhmc/crhmc_walk.hpp @@ -164,6 +164,9 @@ struct CRHMCWalk { inline void disable_adaptive(){ update_modules=false; } + NT get_current_eta() const { + return solver->get_eta(); + } inline void apply(RandomNumberGenerator &rng, int walk_length = 1, bool metropolis_filter = true) diff --git a/include/random_walks/gaussian_helpers.hpp b/include/random_walks/gaussian_helpers.hpp index d1ab5429a..76a151130 100644 --- a/include/random_walks/gaussian_helpers.hpp +++ b/include/random_walks/gaussian_helpers.hpp @@ -9,7 +9,14 @@ NT eval_exp(Point const& p, NT const& a) { return std::exp(-a * p.squared_length()); } - +template +NT eval_exp(Point const& p, MT const& inv_covariance_matrix, NT const& a_next, NT const& a_curr) +{ + Eigen::Matrix dist_vector = p.getCoefficients(); + NT mahalanobis_dist = dist_vector.transpose() * inv_covariance_matrix * dist_vector; + NT log_ratio = (a_curr - a_next) * mahalanobis_dist; + return std::exp(log_ratio); +} template NT get_max(Point const& l, Point const& u, NT const& a_i) diff --git a/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp b/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp new file mode 100644 index 000000000..fd5ce185d --- /dev/null +++ b/include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp @@ -0,0 +1,594 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2024 Vissarion Fisikopoulos +// Copyright (c) 2018-2024 Apostolos Chalkis +// Copyright (c) 2024 Vladimir Necula + +// Contributed and/or modified by Vladimir Necula, as part of Google Summer of +// Code 2024 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +#ifndef VOLUME_COOLING_NON_SPHERICAL_GAUSSIANS_CRHMC_HPP +#define VOLUME_COOLING_NON_SPHERICAL_GAUSSIANS_CRHMC_HPP + +//#define VOLESTI_DEBUG + +#include "volume/volume_cooling_gaussians.hpp" +#include "preprocess/crhmc/crhmc_problem.h" +#include "preprocess/max_inscribed_ellipsoid.hpp" +#include "preprocess/inscribed_ellipsoid_rounding.hpp" + +////////////////////////////// Algorithms + +// Gaussian Anealling + +// Compute a_{i+1} when a_i is given + +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); + +} + + +template +< + int simdLen, + typename Polytope, + typename Point, + typename NT, + typename MT, + 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, + 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_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); + + 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; + const NT tol = 0.00001; + bool done = false; + std::vector fn(N, NT(0.0)); + std::list randPoints; + + // sample N points + PushBackWalkPolicy push_back_policy; + bool raw_output = false; + typedef CrhmcRandomPointGenerator CRHMCRandomPointGenerator; + + CRHMCRandomPointGenerator::apply(problem, p, N, walk_length, randPoints, + push_back_policy, rng, g, f, params, 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, inv_covariance_matrix, new_a, 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 the new a value as a scalar + return last_a * std::pow(ratio, k); +} + +// Compute the sequence of non spherical gaussians +template< + int simdLen, + typename Point, + typename Polytope, + typename NT, + typename MT, + typename RandomNumberGenerator +> +void compute_annealing_schedule(Polytope Pin_copy, + Polytope P_copy, + NT const& ratio, + NT const& C, + NT const& frac, + unsigned int const& N, + unsigned int const& walk_length, + NT const& error, + std::vector& a_vals, + MT const& inv_covariance_matrix, + RandomNumberGenerator& rng) +{ + + 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_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; + + // compute the first variance + // for this we need P + 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 + + //for the rest we need Pin + NT a_stop = 0.0; + const NT tol = 0.001; + unsigned int it = 0; + unsigned int n = Pin_copy.dimension(); + const unsigned int totalSteps = ((int)150/((1.0 - frac) * error)) + 1; + + if (a_vals[0] initial_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); + initial_problem.options.simdLen = simdLen; + crhmc_walk_params initial_params(initial_input.df, initial_p.dimension(), initial_problem.options); + CRHMCWalkType initial_walk = CRHMCWalkType(initial_problem, initial_p, initial_input.df, initial_input.f, initial_params); + + burnIn( + initial_p, walk_length, rng, initial_g, initial_f, initial_params, initial_problem, initial_walk, 2.5 * N); + + //fix eta + initial_eta = initial_walk.get_current_eta(); + + //fix point + start_point = initial_problem.center; + + 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; + + //TODO: potential problem creation and preprocessing optimization + + // Create the CRHMC problem for this variance + int dimension = Pin_copy.dimension(); + func_params f_params = func_params(Point(dimension), a_vals[it], -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + ZeroFunctor zerof; + + Input input = convert2crhmc_input>(Pin_copy, f, g, zerof); + + typedef crhmc_problem CrhmcProblem; + + CrhmcProblem problem = CrhmcProblem(input); + + if(problem.terminate) { return; } + + problem.options.simdLen = simdLen; + + crhmc_walk_params params(input.df, start_point.dimension(), problem.options); + + //eta initialization with the previous walk eta + problem.options.etaInitialize = false; + params.eta = initial_eta; + + // Create the walk object for this problem + CRHMCWalkType walk = CRHMCWalkType(problem, start_point, input.df, input.f, params); + + // Compute some ratios to decide if this is the last gaussian + for (unsigned int j = 0; j < steps; j++) + { + walk.apply(rng, walk_length); + start_point = walk.getPoint(); + curr_its += 1.0; + curr_fn += eval_exp(start_point, inv_covariance_matrix, next_a, a_vals[it]); + } + + //restore the new eta and start point, by looking at the walk after its operations + initial_eta = walk.get_current_eta(); + + // Remove the last gaussian. + // Set the last a_i equal to zero + if (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; + } + } +} + + +template +struct non_gaussian_annealing_parameters +{ + non_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))) + , W(6 * d * d + 800) + {} + + NT frac; + NT ratio; + NT C; + unsigned int N; + unsigned int W; +}; + + +template +< + typename Polytope, + typename RandomNumberGenerator, + typename WalkTypePolicy = CRHMCWalk, + int simdLen = 8 +> +double non_spherical_crhmc_volume_cooling_gaussians(Polytope& Pin, + RandomNumberGenerator& rng, + double const& error = 0.1, + unsigned int const& walk_length = 1) +{ + typedef typename Polytope::PointType Point; + typedef typename Point::FT NT; + typedef typename Polytope::VT VT; + typedef typename Polytope::MT MT; + typedef typename NonSphericalGaussianFunctor::FunctionFunctor Func; + typedef typename NonSphericalGaussianFunctor::GradientFunctor Grad; + typedef typename NonSphericalGaussianFunctor::parameters func_params; + + typedef crhmc_input> Input; + typedef crhmc_problem CrhmcProblem; + + typedef ImplicitMidpointODESolver Solver; + + typedef typename WalkTypePolicy::template Walk + < + Point, + CrhmcProblem, + RandomNumberGenerator, + Grad, + Func, + Solver + > CRHMCWalkType; + + typedef typename WalkTypePolicy::template parameters + < + NT, + Grad + > crhmc_walk_params; + + typedef HPolytope HPOLYTOPE; + + HPOLYTOPE P(Pin.dimension(), Pin.get_mat(), Pin.get_vec()); + HPOLYTOPE newPin(Pin.dimension(), Pin.get_mat(), Pin.get_vec()); + unsigned int n = P.dimension(); + unsigned int m = P.num_of_hyperplanes(); + + //compute inscribed ellipsoid + NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0); + unsigned int maxiter = 100; + P.normalize(); + VT x0 = compute_feasible_point(P.get_mat(), P.get_vec()); + auto ellipsoid_result = compute_inscribed_ellipsoid(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg); + + // extract the covariance matrix and the center of the ellipsoid + MT inv_covariance_matrix = std::get<0>(ellipsoid_result); //this is the covariance to use in the telescopic product + VT center = std::get<1>(ellipsoid_result); + MT covariance_matrix = inv_covariance_matrix.inverse(); + + newPin.shift(center); //we shift the initial polytope so that the origin is the center of the gaussians + P.shift(center); + + // we apply the rounding transformation + Eigen::LLT lltOfA(covariance_matrix); + auto L = lltOfA.matrixL(); + P.linear_transformIt(L); + + + // Initialize the gaussian_annealing_parameters struct + non_gaussian_annealing_parameters parameters(P.dimension()); + + // 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(newPin, P, ratio, C, parameters.frac, N, walk_length, error, a_vals, inv_covariance_matrix, 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); + + MT scaled_matrix = a_vals[0] * inv_covariance_matrix; + NT det_scaled = scaled_matrix.determinant(); + NT vol = std::pow(M_PI, n / 2.0) / std::sqrt(det_scaled); + + unsigned int i=0; + + typedef typename std::vector::iterator viterator; + viterator itsIt = its.begin(); + auto 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 = newPin.dimension(); + func_params f_params = func_params(Point(dimension), *avalsIt, -1, inv_covariance_matrix); + + Func f(f_params); + Grad g(f_params); + + ZeroFunctor zerof; + + //create the crhmc problem + Input input = convert2crhmc_input>(newPin, f, g, zerof); + + typedef crhmc_problem CrhmcProblem; + CrhmcProblem problem = CrhmcProblem(input); + Point p = problem.center; + + if(problem.terminate){return 0;} + problem.options.simdLen=simdLen; + + //create the walk and do the burnIn + crhmc_walk_params params(input.df, p.dimension(), problem.options); + CRHMCWalkType walk = CRHMCWalkType(problem, p, input.df, input.f, params); + + burnIn( + p, walk_length, rng, g, f, params, problem, walk, 2.5 * N); + + while (!done || (*itsIt) < min_steps) + { + walk.apply(rng, walk_length); + p = walk.getPoint(); + *itsIt = *itsIt + 1.0; + + *fnIt += eval_exp(p, inv_covariance_matrix, *(avalsIt+1), *avalsIt); + + NT val = (*fnIt) / (*itsIt); + + last_W[index] = val; + if (val <= min_val) + { + min_val = val; + min_index = index; + } else if (min_index == index) + { + minmaxIt = std::min_element(last_W.begin(), last_W.end()); + min_val = *minmaxIt; + min_index = std::distance(last_W.begin(), minmaxIt); + } + + if (val >= 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 = "<