Rename LinearProgramFunctor to ExponentialFunctor, add inverse varian…

…ce parameter, add define lp_solve guards to orderpolytope
GeomScale · Feb 14, 2024 · f2ba8f5 · f2ba8f5
1 parent e767f46
commit f2ba8f5
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 5 deletions.
diff --git 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
@@ -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);
     }
 
   };