diff --git a/examples/sample_points/simple_example.cpp b/examples/sample_points/simple_example.cpp index 5f11acf2c..25f7e1e91 100644 --- a/examples/sample_points/simple_example.cpp +++ b/examples/sample_points/simple_example.cpp @@ -25,12 +25,12 @@ using NT = double; using MT = Eigen::Matrix; using VT = Eigen::Matrix; -template -void sample_points_eigen_matrix(HPolytopeType const& HP, Point const& q, Walk const& walk, +template +void sample_points_eigen_matrix(PolytopeOrProblem const& HP, Point const& q, Walk const& walk, Distribution const& distr, RNGType rng, int walk_len, int rnum, int nburns) { - MT samples(HP.dimension(), rnum); + MT samples(q.dimension(), rnum); sample_points(HP, q, walk, distr, rng, walk_len, rnum, nburns, samples); @@ -96,84 +96,48 @@ struct CustomFunctor { }; }; -struct CustomGaussianFunctor { - - template - < - typename NT, - typename Point - > - struct parameters { - Point x0; - NT a; - NT eta; - unsigned int order; - NT L; // Lipschitz constant for gradient - NT m; // Strong convexity constant - NT kappa; // Condition number - - parameters(Point x0_, NT a_, NT eta_) : - x0(x0_), a(a_), eta(eta_), order(2), L(2 * a_), m(2 * a_), kappa(1) {}; - - }; - - template - struct GradientFunctor { - typedef typename Point::FT NT; - typedef std::vector pts; - - parameters ¶ms; - - GradientFunctor(parameters ¶ms_) : params(params_) {}; - - // The index i represents the state vector index - /* - Point operator() (unsigned int const& i, pts const& xs, NT const& t) const { - if (i == params.order - 1) { - Point y = (-2.0 * params.a) * (xs[0] - params.x0); - return y; - } else { - return xs[i + 1]; // returns derivative - } - }*/ - Point operator()(Point const&x){ - Point y = (-2.0 * params.a) * (x - params.x0); - return y; - } - }; - - template - struct FunctionFunctor { - typedef typename Point::FT NT; - - parameters ¶ms; - - FunctionFunctor(parameters ¶ms_) : params(params_) {}; - - // The index i represents the state vector index - NT operator() (Point const& x) const { - Point y = x - params.x0; - return params.a * y.dot(y); - } - }; - - template - struct HessianFunctor { - typedef typename Point::FT NT; +inline bool exists_check(const std::string &name) { + std::ifstream f(name.c_str()); + return f.good(); +} - parameters ¶ms; +template< typename SpMat, typename VT> +void load_crhmc_problem(SpMat &A, VT &b, VT &lb, VT &ub, int &dimension, + std::string problem_name) { + { + std::string fileName("../crhmc_sampling/data/"); + fileName.append(problem_name); + fileName.append(".mm"); + if(!exists_check(fileName)){ + std::cerr<<"Problem does not exist.\n"; + exit(1);} + SpMat X; + loadMarket(X, fileName); + int m = X.rows(); + dimension = X.cols() - 1; + A = X.leftCols(dimension); + b = VT(X.col(dimension)); + } + { + std::string fileName("../crhmc_sampling/data/"); + fileName.append(problem_name); + fileName.append("_bounds.mm"); + if(!exists_check(fileName)){ + std::cerr<<"Problem does not exist.\n"; + exit(1);} + SpMat bounds; + loadMarket(bounds, fileName); + lb = VT(bounds.col(0)); + ub = VT(bounds.col(1)); + } +} - HessianFunctor(parameters ¶ms_) : params(params_) {}; - // The index i represents the state vector index - Point operator() (Point const& x) const { - return (2.0 * params.a) * Point::all_ones(x.dimension()); - } - }; +int main() { + // NEW INTERFACE Sampling -}; + // Inputs: -int main() { // Generating a 3-dimensional cube centered at origin HPolytopeType HP = generate_cube(10, false); std::cout<<"Polytope: \n"; @@ -185,7 +149,19 @@ int main() { Point q(HP.dimension()); RNGType rng(HP.dimension()); - // NEW INTERFACE Sampling + // Generating a sparse polytope/problem + using SpMat = Eigen::SparseMatrix; + using ConstraintProblem =constraint_problem; + std::string problem_name("simplex3"); + std::cerr << "CRHMC on " << problem_name << "\n"; + SpMat As; + VT b, lb, ub; + int dimension; + load_crhmc_problem(As, b, lb, ub, dimension, problem_name); + ConstraintProblem problem = ConstraintProblem(dimension); + problem.set_equality_constraints(As, b); + problem.set_bounds(lb, ub); + // Walks AcceleratedBilliardWalk abill_walk; @@ -245,15 +221,16 @@ int main() { NegativeLogprobFunctorR fr(params_r); LogConcaveDistribution logconcave_reflective(gr, fr, params_r.L); - using NegativeGradientFunctor = CustomGaussianFunctor::GradientFunctor; - using NegativeLogprobFunctor = CustomGaussianFunctor::FunctionFunctor; - using HessianFunctor = CustomGaussianFunctor::HessianFunctor; - CustomGaussianFunctor::parameters params(x0, 0.5, 1); + using NegativeGradientFunctor = GaussianFunctor::GradientFunctor; + using NegativeLogprobFunctor = GaussianFunctor::FunctionFunctor; + using HessianFunctor = GaussianFunctor::HessianFunctor; + GaussianFunctor::parameters params(x0, 0.5, 1); NegativeGradientFunctor g(params); NegativeLogprobFunctor f(params); HessianFunctor h(params); LogConcaveDistribution logconcave_crhmc(g, f, h, params.L); + LogConcaveDistribution logconcave_ref_gaus(g, f, params.L); // Sampling @@ -293,8 +270,10 @@ int main() { std::cout << "logconcave" << std::endl; sample_points_eigen_matrix(HP, q, hmc_walk, logconcave_reflective, rng, walk_len, rnum, nburns); sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_reflective, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(HP, q, nhmc_walk, logconcave_ref_gaus, rng, walk_len, rnum, nburns); sample_points_eigen_matrix(HP, q, crhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns); + sample_points_eigen_matrix(problem, q, crhmc_walk, logconcave_crhmc, rng, walk_len, rnum, nburns); std::cout << "fix the following" << std::endl; diff --git a/include/sampling/sample_points.hpp b/include/sampling/sample_points.hpp index 0cf121bc4..48e839b9a 100644 --- a/include/sampling/sample_points.hpp +++ b/include/sampling/sample_points.hpp @@ -258,9 +258,8 @@ void sample_points(Polytope& P, // TODO: make it a const& using HPolytope = typename std::remove_const::type; HPolytope HP = P; //TODO: avoid the copy - constexpr int simdLen = 8; + constexpr int simdLen = 8; //TODO: input parameter using NT = double; - using MT = typename HPolytope::MT; int dimension = HP.dimension(); @@ -298,19 +297,19 @@ void sample_points(Polytope& P, // TODO: make it a const& >; using Walk = typename WalkType::template Walk - < - Point, - CrhmcProblem, - RandomNumberGenerator, - NegativeGradientFunctor, - NegativeLogprobFunctor, - Solver - >; + < + Point, + CrhmcProblem, + RandomNumberGenerator, + NegativeGradientFunctor, + NegativeLogprobFunctor, + Solver + >; using WalkParams = typename WalkType::template parameters - < - NT, - NegativeGradientFunctor - >; + < + NT, + NegativeGradientFunctor + >; Point p = Point(problem.center); problem.options.simdLen=simdLen; WalkParams params(distribution.L, p.dimension(), problem.options); @@ -333,7 +332,7 @@ void sample_points(Polytope& P, // TODO: make it a const& { walk.apply(rng, walk_len); if (walk.P.terminate) {return;} - MT x = raw_output ? walk.x : walk.getPoints(); + auto x = raw_output ? walk.x : walk.getPoints(); if ((i + 1) * simdLen > rnum) {