Fix 1913; initialization bug; dimension bug; templates (#1914)

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/linprog.cpp

@@ -10,7 +10,6 @@
 #include "slice_into.h"
 #include "find.h"
 #include "colon.h"
-#include <iostream>
 
 //#define IGL_LINPROG_VERBOSE
 IGL_INLINE bool igl::linprog(
@@ -242,7 +241,7 @@ IGL_INLINE bool igl::linprog(
     s*=1e6*c.array().abs().maxCoeff();
     s.head(n) = c;
   }
-  x.resize(std::max(B.maxCoeff()+1,n));
+  x.setZero(std::max(B.maxCoeff()+1,n));
   igl::slice_into(xb,B,x);
   x = x.head(n).eval();
   return success;
@@ -281,9 +280,12 @@ IGL_INLINE bool igl::linprog(
   MatrixXd BSP(0,2*n+m);
   if(p>0)
   {
-    MatrixXd BS(p,2*n);
-    BS<<B,MatrixXd::Zero(p,n);
+    // B ∈ ℝ^(p × n)
+    MatrixXd BS(p,n+m);
+    BS<<B,MatrixXd::Zero(p,m);
+    // BS ∈ ℝ^(p × n+m)
     BSP = BS*P;
+    // BSP ∈ ℝ^(p × 2n+m)
   }
 
   VectorXd fSP = VectorXd::Ones(2*n+m);
@@ -296,7 +298,24 @@ IGL_INLINE bool igl::linprog(
   bb<<bS,c;
 
   VectorXd xxs;
+  // min   ccᵀxxs
+  // s.t.  AA xxs ≤ bb
+  //          xxs ≥ 0
+  //        
+  // x = x⁺ - x⁻
+  //
+  //    P
+  // .--^---.
+  // [I -I 0  [x⁺   = [x
+  //  0  0 I]  x⁻      s]
+  //           s]
+  // Pᵀ [xᵀ sᵀ] = xxsᵀ
+  //
+  // min  [fᵀ -fᵀ 𝟙ᵀ] [x⁺;x⁻;s]
+  // s.t.  AA [x⁺;x⁻;s] ≤ b
+  // s.t.  [x⁺;x⁻;s] ≥ 0
   bool ret = linprog(cc,AA,bb,0,xxs);
+  // x = P(1:n,:) xxs
   x = P.block(0,0,n,2*n+m)*xxs;
   return ret;
 }

include/igl/linprog.h

@@ -13,7 +13,7 @@ namespace igl
 {
   // Solve a linear program given in "standard form"
   //
-  // min  f'x
+  // min  c'x
   // s.t. A(    1:k,:) x <= b(1:k)
   //      A(k+1:end,:) x = b(k+1:end)
   //   ** x >= 0 **
@@ -33,7 +33,7 @@ namespace igl
     const Eigen::MatrixXd & A,
     const Eigen::VectorXd & b,
     const int k,
-    Eigen::VectorXd & f);
+    Eigen::VectorXd & x);
   
   // Wrapper in friendlier general form (no implicit bounds on x)
   //

include/igl/slice.cpp

@@ -188,6 +188,8 @@ IGL_INLINE DerivedX igl::slice(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Matrix<double, -1, 1, 0, -1, 1> igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&);
+template Eigen::Matrix<int, -1, 1, 0, -1, 1> igl::slice<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::DenseBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1>, -1, 1, true> > const&, int, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
 // generated by autoexplicit.sh

tests/include/igl/linprog.cpp

@@ -0,0 +1,41 @@
+#include <test_common.h>
+#include <igl/linprog.h>
+
+TEST_CASE("linprog: 2D-inequality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(0,2);
+  Eigen::VectorXd c = Eigen::VectorXd(0);
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<6,4;
+  test_common::assert_near( x, x_correct, 1e-10);
+}
+
+TEST_CASE("linprog: 2D-inequality+2-equality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(2,2);B<<1,0,0,1;
+  Eigen::VectorXd c = Eigen::VectorXd(2);c<<4,4;
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<4,4;
+  test_common::assert_near( x, x_correct, 1e-10);
+}
+
+TEST_CASE("linprog: 2D-inequality+1-equality", "[igl]" )
+{
+  Eigen::VectorXd f(2); f << -2, -1;
+  Eigen::VectorXd b(3); b << 8, 14, 4;
+  Eigen::MatrixXd A(3, 2); A << 2, -1, 1, 2, -1, 1;
+  Eigen::MatrixXd B = Eigen::MatrixXd(1,2);B<<1,0;
+  Eigen::VectorXd c = Eigen::VectorXd(1);c<<4;
+  Eigen::VectorXd x;
+  igl::linprog(f, A, b, B, c, x);
+  Eigen::VectorXd x_correct(2); x_correct<<4,5;
+  test_common::assert_near( x, x_correct, 1e-10);
+}