diff --git a/R-proj/src/Makevars b/R-proj/src/Makevars index b642e2fbe..ea3607acf 100644 --- a/R-proj/src/Makevars +++ b/R-proj/src/Makevars @@ -2,8 +2,6 @@ PKG_CPPFLAGS= -I../../external/boost -I../../external/LPsolve_src/run_headers -I PKG_CXXFLAGS= -Wno-deprecated-declarations -lm -ldl -Wno-ignored-attributes -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES -CXX_STD = CXX11 - PKG_LIBS=-LRproj_externals/lp_solve -llp_solve -L../../external/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) $(SHLIB): Rproj_externals/lp_solve/liblp_solve.a ../../external/PackedCSparse/qd/libqd.a diff --git a/R-proj/src/Makevars.win b/R-proj/src/Makevars.win index 17d28e206..83ed4e0c6 100644 --- a/R-proj/src/Makevars.win +++ b/R-proj/src/Makevars.win @@ -1,6 +1,5 @@ PKG_CPPFLAGS=-I../../external/boost -I../../external/LPsolve_src/run_headers -I../../external/minimum_ellipsoid -I../../include -I../../include/convex_bodies/spectrahedra PKG_CXXFLAGS= -Wno-deprecated-declarations -lm -ldl -Wno-ignored-attributes -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES -CXX_STD = CXX11 PKG_LIBS=-LRproj_externals/lp_solve -llp_solve -L../../external/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/cran_gen/Makevars b/cran_gen/Makevars index 28c42dcad..8194d4c8a 100644 --- a/cran_gen/Makevars +++ b/cran_gen/Makevars @@ -1,6 +1,5 @@ PKG_CPPFLAGS=-Iexternal -Iexternal/lpsolve/headers/run_headers -Iexternal/minimum_ellipsoid -Iinclude -Iinclude/convex_bodies/spectrahedra PKG_CXXFLAGS= -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES -CXX_STD = CXX11 PKG_LIBS=-Lexternal/lpsolve/build/lp_solve -llp_solve -Lexternal/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/cran_gen/Makevars.win b/cran_gen/Makevars.win index b93299131..e4f79d991 100644 --- a/cran_gen/Makevars.win +++ b/cran_gen/Makevars.win @@ -1,6 +1,5 @@ PKG_CPPFLAGS=-Iexternal -Iexternal/lpsolve/headers/run_headers -Iexternal/minimum_ellipsoid -Iinclude -Iinclude/convex_bodies/spectrahedra PKG_CXXFLAGS= -lm -ldl -DBOOST_NO_AUTO_PTR -DDISABLE_NLP_ORACLES -CXX_STD = CXX11 PKG_LIBS=-Lexternal/lpsolve/build/lp_solve -llp_solve -Lexternal/PackedCSparse/qd -lqd $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) diff --git a/include/convex_bodies/orderpolytope.h b/include/convex_bodies/orderpolytope.h index 66854234b..2df5c6619 100644 --- a/include/convex_bodies/orderpolytope.h +++ b/include/convex_bodies/orderpolytope.h @@ -208,7 +208,32 @@ class OrderPolytope { std::pair ComputeInnerBall() { normalize(); - return ComputeChebychevBall(_A, b); + std::pair inner_ball; + #ifndef DISABLE_LPSOLVE + inner_ball = ComputeChebychevBall(_A, b); // use lpsolve library + #else + + if (inner_ball.second <= NT(0)) { + + NT const tol = 0.00000001; + std::tuple inner_ball = max_inscribed_ball(_A, b, 150, tol); + + // check if the solution is feasible + if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < NT(0) || + std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) || + !std::get<2>(inner_ball) || is_inner_point_nan_inf(std::get<0>(inner_ball))) + { + inner_ball.second = -1.0; + } else + { + inner_ball.first = Point(std::get<0>(inner_ball)); + inner_ball.second = std::get<1>(inner_ball); + } + } + #endif + + return inner_ball; + } diff --git a/include/ode_solvers/oracle_functors.hpp b/include/ode_solvers/oracle_functors.hpp index 3aea450b2..35fc8887b 100644 --- a/include/ode_solvers/oracle_functors.hpp +++ b/include/ode_solvers/oracle_functors.hpp @@ -243,7 +243,7 @@ struct IsotropicLinearFunctor { }; -struct LinearProgramFunctor { +struct ExponentialFunctor { // Sample from linear program c^T x (exponential density) template < @@ -256,8 +256,10 @@ struct LinearProgramFunctor { NT m; // Strong convexity constant NT kappa; // Condition number Point c; // Coefficients of LP objective + NT a; // Inverse variance - parameters(Point c_) : order(2), L(1), m(1), kappa(1), c(c_) {}; + parameters(Point c_) : order(2), L(1), m(1), kappa(1), c(c_), a(1.0) {}; + parameters(Point c_, NT a_) : order(2), L(1), m(1), kappa(1), c(c_), a(a_) {}; }; @@ -277,7 +279,7 @@ struct LinearProgramFunctor { Point operator() (unsigned int const& i, pts const& xs, NT const& t) const { if (i == params.order - 1) { Point y(params.c); - return (-1.0) * y; + return (-params.a) * y; } else { return xs[i + 1]; // returns derivative } @@ -298,7 +300,7 @@ struct LinearProgramFunctor { // The index i represents the state vector index NT operator() (Point const& x) const { - return x.dot(params.c); + return params.a * x.dot(params.c); } }; diff --git a/test/logconcave_sampling_test.cpp b/test/logconcave_sampling_test.cpp index 8e53f489b..1f67b2851 100644 --- a/test/logconcave_sampling_test.cpp +++ b/test/logconcave_sampling_test.cpp @@ -334,12 +334,12 @@ void benchmark_nuts_hmc(bool truncated) { bool automatic_burnin = false; std::chrono::time_point start, stop; - for (unsigned int dim = dim_min; dim <= dim_max; dim+=10) + for (unsigned int dim = dim_min; dim <= dim_max; dim+=10) { MT samples(dim, n_samples); Point x0(dim); NutsHamiltonianMonteCarloWalk::parameters hmc_params(F, dim); - if (truncated) + if (truncated) { Hpolytope P = generate_cube(dim, false); @@ -784,8 +784,8 @@ void benchmark_polytope_linear_program_optimization( typedef std::vector pts; typedef boost::mt19937 RNGType; typedef BoostRandomNumberGenerator RandomNumberGenerator; - typedef LinearProgramFunctor::GradientFunctor NegativeGradientFunctor; - typedef LinearProgramFunctor::FunctionFunctor NegativeLogprobFunctor; + typedef ExponentialFunctor::GradientFunctor NegativeGradientFunctor; + typedef ExponentialFunctor::FunctionFunctor NegativeLogprobFunctor; typedef OptimizationFunctor::GradientFunctor NegativeGradientOptimizationFunctor; typedef OptimizationFunctor::FunctionFunctor lp_params(coeffs); + ExponentialFunctor::parameters lp_params(coeffs); NegativeGradientFunctor F_lp(lp_params); NegativeLogprobFunctor f_lp(lp_params); @@ -1116,8 +1116,8 @@ void call_test_exp_sampling() { typedef HPolytope Hpolytope; std::string name; std::vector> polytopes; - - + + if (exists_check("metabolic_full_dim/e_coli_biomass_function.txt") && exists_check("metabolic_full_dim/polytope_e_coli.ine")){ Point biomass_function_e_coli = load_biomass_function("metabolic_full_dim/e_coli_biomass_function.txt"); polytopes.push_back(std::make_tuple(read_polytope("metabolic_full_dim/polytope_e_coli.ine"), biomass_function_e_coli, "e_coli", true));