From 1ee9073569eda1bf31d2aa4752a85ea260ff5151 Mon Sep 17 00:00:00 2001
From: vfisikop <vissarion.fisikopoulos@gmail.com>
Date: Wed, 14 Feb 2024 11:54:00 +0200
Subject: [PATCH] Rename LinearProgramFunctor to ExponentialFunctor, add
 inverse variance parameter, add define lp_solve guards to orderpolytope

---
 include/convex_bodies/orderpolytope.h   | 27 ++++++++++++++++++++++++-
 include/ode_solvers/oracle_functors.hpp | 10 +++++----
 test/logconcave_sampling_test.cpp       | 14 ++++++-------
 3 files changed, 39 insertions(+), 12 deletions(-)

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<Point, NT> ComputeInnerBall()
     {
         normalize();
-        return ComputeChebychevBall<NT, Point>(_A, b);
+        std::pair<Point, NT> inner_ball;
+        #ifndef DISABLE_LPSOLVE
+            inner_ball = ComputeChebychevBall<NT, Point>(_A, b); // use lpsolve library
+        #else
+
+            if (inner_ball.second <= NT(0)) {
+
+                NT const tol = 0.00000001;
+                std::tuple<VT, NT, bool> 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<std::chrono::high_resolution_clock> 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<NT, NegativeGradientFunctor> hmc_params(F, dim);
-      if (truncated) 
+      if (truncated)
       {
         Hpolytope P = generate_cube<Hpolytope>(dim, false);
 
@@ -784,8 +784,8 @@ void benchmark_polytope_linear_program_optimization(
     typedef std::vector<Point> pts;
     typedef boost::mt19937 RNGType;
     typedef BoostRandomNumberGenerator<RNGType, NT> RandomNumberGenerator;
-    typedef LinearProgramFunctor::GradientFunctor<Point> NegativeGradientFunctor;
-    typedef LinearProgramFunctor::FunctionFunctor<Point> NegativeLogprobFunctor;
+    typedef ExponentialFunctor::GradientFunctor<Point> NegativeGradientFunctor;
+    typedef ExponentialFunctor::FunctionFunctor<Point> NegativeLogprobFunctor;
     typedef OptimizationFunctor::GradientFunctor<Point, NegativeLogprobFunctor,
         NegativeGradientFunctor> NegativeGradientOptimizationFunctor;
     typedef OptimizationFunctor::FunctionFunctor<Point, NegativeLogprobFunctor,
@@ -816,7 +816,7 @@ void benchmark_polytope_linear_program_optimization(
     }
 
     // Declare oracles for LP
-    LinearProgramFunctor::parameters<NT, Point> lp_params(coeffs);
+    ExponentialFunctor::parameters<NT, Point> lp_params(coeffs);
 
     NegativeGradientFunctor F_lp(lp_params);
     NegativeLogprobFunctor f_lp(lp_params);
@@ -1116,8 +1116,8 @@ void call_test_exp_sampling() {
     typedef HPolytope<Point> Hpolytope;
     std::string name;
     std::vector<std::tuple<Hpolytope, Point, std::string, bool>> 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<Point, NT>("metabolic_full_dim/e_coli_biomass_function.txt");
       polytopes.push_back(std::make_tuple(read_polytope<Hpolytope, NT>("metabolic_full_dim/polytope_e_coli.ine"), biomass_function_e_coli, "e_coli", true));