libigl/tests/include/igl/quadprog.cpp

#include <test_common.h>
#include <igl/quadprog.h>
#include <igl/EPS.h>


TEST_CASE("quadprog: linear inequalities", "[igl]" )
{
  const Eigen::Matrix<double,6,6> H = (Eigen::Matrix<double,6,6>(6,6)<<111,126,120,107,165,130,126,172,134,148,192,163,120,134,164,131,177,144,107,148,131,207,205,169,165,192,177,205,310,219,130,163,144,169,219,190).finished();
  const Eigen::Matrix<double,6,1> f = (Eigen::Matrix<double,6,1>(6,1)<<3,1,1,9,6,5).finished();
  const Eigen::Matrix<double,2,6> A = (Eigen::Matrix<double,2,6>(2,6)<<8,4,2,0,0,8,7,7,7,3,1,7).finished();
  const Eigen::Matrix<double,2,1> lbi = (Eigen::Matrix<double,2,1>(2,1)<<-0.5,-6).finished();
  const Eigen::Matrix<double,2,1> ubi = (Eigen::Matrix<double,2,1>(2,1)<<0.5,6).finished();
  const Eigen::Matrix<double,6,1> lb = (Eigen::Matrix<double,6,1>(6,1)<<-0.1,-0.1,-0.1,-0.1,-0.1,-0.1).finished();
  const Eigen::Matrix<double,6,1> ub = (Eigen::Matrix<double,6,1>(6,1)<<0.1,0.1,0.1,0.1,0.1,0.1).finished();
  Eigen::Matrix<double,6,1> x = igl::quadprog(H,f,A,lbi,ubi,lb,ub);
  //std::cout<<igl::matlab_format(H,"H")<<std::endl;
  //std::cout<<igl::matlab_format(f,"f")<<std::endl;
  //std::cout<<igl::matlab_format(A,"A")<<std::endl;
  //std::cout<<igl::matlab_format(lbi,"lbi")<<std::endl;
  //std::cout<<igl::matlab_format(ubi,"ubi")<<std::endl;
  //std::cout<<igl::matlab_format(lb,"lb")<<std::endl;
  //std::cout<<igl::matlab_format(ub,"ub")<<std::endl;
  //std::cout<<igl::matlab_format(x,"x")<<std::endl;
  REQUIRE(abs(x(0)- -0.1)<1e-10);
  REQUIRE(abs(x(1)-  0.1)<1e-10);
  REQUIRE(abs(x(2)-  0.07227715420917860)<1e-10);
  REQUIRE(abs(x(3)- -0.1)<1e-10);
  REQUIRE(abs(x(4)-  0.01839231579976747)<1e-10);
  REQUIRE(abs(x(5)- -0.03056928855229465)<1e-10);
}

TEST_CASE("quadprog: box3", "[igl]" )
{
  {
    Eigen::Matrix3d H = (Eigen::Matrix3d(3,3)<<
        0.240548455386281,   0.237314308102107,  0.0436993831501944,
        0.237314308102107,   0.326254049041854,-0.00021896091952631,
        0.0436993831501944,-0.00021896091952631,   0.171175681280756
        ).finished();
    Eigen::Vector3d f = (Eigen::Vector3d(3,1)<<
        0.222182612270718,
        0.503254616893693,
        -0.184619987497072
        ).finished();
    Eigen::Vector3d lb = (Eigen::Vector3d(3,1)<<
        -1,
        -1,
        -1
        ).finished();
    Eigen::Vector3d ub = (Eigen::Vector3d(3,1)<<
        1,
        1,
        1
        ).finished(); 
    Eigen::Vector3d x = igl::quadprog(H,f,lb,ub);
    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
    REQUIRE(abs(x(0)- -0.118760635036839)<1e-7);
    REQUIRE(abs(x(1)- -1)<1e-7);
    REQUIRE(abs(x(2)- +1)<1e-7);
  }
  {
    Eigen::Matrix3d H;
    H<<1,0,0,0,1,0,0,0,1;
    Eigen::Vector3d f( 0.5,-0.5,-0.5);
    Eigen::Vector3d lb(0,0,0);
    Eigen::Vector3d ub(1,1,1);
    Eigen::Vector3d x = igl::quadprog(H,f,lb,ub);
    REQUIRE(x(0)==0.0);
    REQUIRE(x(1)==0.5);
    REQUIRE(x(1)==0.5);
  }
  {
    Eigen::Matrix3d H = (Eigen::Matrix3d(3,3)<<
        1.06020279605748, 0.387347953430924,-0.653587847224834,
        0.387347953430924, 0.323001970642631, -0.80259721932688,
        -0.653587847224834, -0.80259721932688,  2.73709523329989
        ).finished();
    Eigen::Vector3d f = (Eigen::Vector3d(3,1)<<
        -0.0155804986732503,
        -0.00161174173383921,
        -0.00903647945917485
        ).finished();
    Eigen::Vector3d lb = (Eigen::Vector3d(3,1)<<
        0,
        0,
        0
        ).finished();
    Eigen::Vector3d ub = (Eigen::Vector3d(3,1)<<
        0.015625,
        0.015625,
        0.015625
        ).finished();
    Eigen::Vector3d x = igl::quadprog(H,f,lb,ub);
    Eigen::Vector3d xexact(0.015625, 0.013732474124087, 0.0110593284260843);
    REQUIRE((x-xexact).array().abs().maxCoeff() < 1e-4);
  }
}

TEST_CASE("quadprog: box2", "[igl]" )
{
  {
    Eigen::Matrix2d H = (Eigen::Matrix2d(2,2)<<
        0.683698654982294,-0.0521997092763332,
        -0.0521997092763332,  0.800535063458999
        ).finished();
    Eigen::Vector2d f = (Eigen::Vector2d(2,1)<<
        -0.733716332234936,
        2.56401312736278
        ).finished();
    Eigen::Vector2d lb = (Eigen::Vector2d(2,1)<<
        -1,
        -1
        ).finished();
    Eigen::Vector2d ub = (Eigen::Vector2d(2,1)<<
        1,
        1
        ).finished();
  
    Eigen::Vector2d x = igl::quadprog(H,f,lb,ub);
    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
    REQUIRE(abs(x(0)-0.99680848864073357)<1e-7);
    REQUIRE(abs(x(1)- -1.)<1e-7);
  }
  {
    Eigen::Matrix2d H = (Eigen::Matrix2d(2,2)<<
        0.0708025678527926,-0.158030756288795,
        -0.158030756288795, 0.360468620664163
        ).finished();
    Eigen::Vector2d f = (Eigen::Vector2d(2,1)<<
        0.207229875768196,
        -0.547595351878845
        ).finished();
    Eigen::Vector2d lb = (Eigen::Vector2d(2,1)<<
        -1,
        -1
        ).finished();
    Eigen::Vector2d ub = (Eigen::Vector2d(2,1)<<
        1,
        1
        ).finished();

    Eigen::Vector2d x = igl::quadprog(H,f,lb,ub);
    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
    REQUIRE(abs(x(0)- -0.69487761491492939)<1e-7);
    REQUIRE(abs(x(1)-1.0)<1e-7);
  }
  {
    Eigen::Matrix2d H;
    H<<0.4000,-0.2000,-0.2000,1.0000;
    Eigen::Vector2d f(-0.3000,4.0000);
    Eigen::Vector2d lb(0,0);
    Eigen::Vector2d ub(1,1);
    Eigen::Vector2d x = igl::quadprog(H,f,lb,ub);
    //std::cout<<igl::matlab_format(x,"x")<<std::endl;
    REQUIRE(abs(x(0)-0.75)<2e-16);
    REQUIRE(abs(x(1)-0.0)<2e-16);
  }
}