#include <CGAL/QP_solver.h>
#include <CGAL/QP_full_exact_pricing.h>
#include <CGAL/QP_partial_exact_pricing.h>
#include <CGAL/QP_full_filtered_pricing.h>
#include <CGAL/QP_partial_filtered_pricing.h>
#include <CGAL/_QP_solver/Double.h>
#include <iostream>

typedef  std::vector<double>  Vector;
typedef  std::vector<Vector>  Matrix;

typedef CGAL::QP_transform_iterator_1<Matrix::const_iterator>  Vector_iterator;
typedef  Vector::const_iterator                                 Entry_iterator;

struct QPSolverTraits {    
    typedef  GMP::Double  ET;
    typedef  Vector_iterator  A_iterator;
    typedef   Entry_iterator  B_iterator;
    typedef   Entry_iterator  C_iterator;
    typedef  Vector_iterator  D_iterator;

    enum Row_type { LESS_EQUAL, EQUAL, GREATER_EQUAL};
    typedef  Row_type*  Row_type_iterator;
    
    typedef  CGAL::Tag_false  Is_linear;
    typedef  CGAL::Tag_false  Is_symmetric;
    typedef  CGAL::Tag_false  Has_full_row_rank;
    typedef  CGAL::Tag_false  Use_perturbation;
};

typedef CGAL::QP_solver<QPSolverTraits>   Solver;

int main( int argc, char** argv)
{

  Matrix  A( 3);
  Vector  b;
  Vector  c;
  Matrix  D( 3);

  // constraint matrix A, column by column
  A[ 0].push_back( -4); A[ 1].push_back(  2); A[ 2].push_back(  0); 
  A[ 0].push_back(  1); A[ 1].push_back(  1); A[ 2].push_back(  1);
      
  // column vector b
  b.push_back( -8); b.push_back(  2);
      
  // linear part of objective function, vector c
  c.push_back(  0); c.push_back(  5); c.push_back(  0);
      
  // quadratic part of objective function, matrix D, row by row
  D[ 0].push_back(  64); D[ 0].push_back( -16); D[ 0].push_back(   0);
  D[ 1].push_back( -16); D[ 1].push_back(   4); D[ 1].push_back(   2);
  D[ 2].push_back(   0); D[ 2].push_back(  -2); D[ 2].push_back(   0); 

  // constraint types
  QPSolverTraits::Row_type*  row_types = new QPSolverTraits::Row_type[ 2];
  row_types[ 0] = QPSolverTraits::GREATER_EQUAL;
  row_types[ 1] = QPSolverTraits::GREATER_EQUAL;
    
  // solve qp with pricing strategy defaulting to full exact pricing strategy
  Solver        qp (3, 2, Vector_iterator(A.begin()), b.begin(), c.begin(),
                    Vector_iterator(D.begin()), row_types);
				     
  // query solution, if qp is optimal				     
  if (qp.status() == Solver::OPTIMAL) {
    Solver::Variable_value_iterator v_it;
    Solver::Variable_value_iterator e_it = qp.variables_value_end(); 
    for (int i = 0, v_it = qp.variables_value_begin(); v_it != e_it; ++v_it, ++i) {
      std::cout << "x[" << i << "]= " << *v_it << std::endl;
    }
  }
  delete[] row_types; 
  return 0;
}