diff --git a/.gitattributes b/.gitattributes index 65120c52fc9..843a2d1280c 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1642,12 +1642,13 @@ QP_solver/examples/QP_solver/first_lp.cpp -text QP_solver/examples/QP_solver/first_nonnegative_lp.cpp -text QP_solver/examples/QP_solver/first_nonnegative_qp.cpp -text QP_solver/examples/QP_solver/first_qp.cpp -text +QP_solver/examples/QP_solver/important_variables.cpp -text QP_solver/examples/QP_solver/integer_qp_solver.cin -text QP_solver/examples/QP_solver/integer_qp_solver.data -text -QP_solver/examples/QP_solver/is_in_convex_hull.h -text -QP_solver/examples/QP_solver/is_in_convex_hull2.h -text QP_solver/examples/QP_solver/rational_qp_solver.cin -text QP_solver/examples/QP_solver/rational_qp_solver.data -text +QP_solver/examples/QP_solver/solve_convex_hull_lp.h -text +QP_solver/examples/QP_solver/solve_convex_hull_lp2.h -text QP_solver/include/CGAL/QP_functions.h -text QP_solver/include/CGAL/QP_models.h -text QP_solver/include/CGAL/QP_solution.h -text diff --git a/QP_solver/doc_tex/QP_solver/main.tex b/QP_solver/doc_tex/QP_solver/main.tex index 8b15dc556c3..2b05b978978 100644 --- a/QP_solver/doc_tex/QP_solver/main.tex +++ b/QP_solver/doc_tex/QP_solver/main.tex @@ -330,14 +330,19 @@ equivalently, if there are $\mu_1,\ldots,\mu_n$ \] The linear program now tests for the existence of nonnegative $\mu_j$ -that satisfy the latter equation. Below is the code; it defines a function -that performs the test, given $p$ and $p_1,\ldots,p_n$ (through an iterator -range), and then calls it with $p_1=(0,0), p_2=(10,0), p_3=(0,10)$ fixed -(they define a triangle) and all integral points $p$ in $[0,10]^2$. We know -that $p$ is in the convex hull of $\{p_1,p_2,p_3\}$ if and only if its two -coordinates sum up to $10$ at most. +that satisfy the latter equation. Below is the code; it defines a +function that solves the linear program, given $p$ and +$p_1,\ldots,p_n$ (through an iterator range). -\ccIncludeExampleCode{QP_solver/is_in_convex_hull.cpp} +\ccIncludeExampleCode{QP_solver/solve_convex_hull_lp.h} + +To see this in action, let us call it with $p_1=(0,0), p_2=(10,0), +p_3=(0,10)$ fixed (they define a triangle) and all integral points $p$ +in $[0,10]^2$. We know that $p$ is in the convex hull of +$\{p_1,p_2,p_3\}$ if and only if its two coordinates sum up to $10$ at +most. + +\ccIncludeExampleCode{QP_solver/convex_hull_containment.cpp} \subsection{Using Makers} You already noticed in the previous example that the actual @@ -354,9 +359,25 @@ You can avoid the explicit construction of the type if you only need an expression of it, e.g.\ to pass it directly as an argument to the solving function. Here is how this works. -\ccIncludeExampleCode{QP_solver/is_in_convex_hull2.cpp} +\ccIncludeExampleCode{QP_solver/solve_convex_hull_lp2.h} \section{The Important Variables} +If we have an optimal solution $x^*$ of a linear or quadratic program, +the ``important'' variables are typically the ones that are not on +their bounds. In case of a nonnegative program, these are the nonzero +variables. Going back to the example of the previous Section +\ref{sec:QP-working_from_iterators}, we can easily interpret their +importance: the nonzero variables correspond to points $p_j$ that +actually contribute to the convex combination that yields $p$. + +The following example shows how we can access the important variables. +We generate a set of points on a line and then find the ones that +contribute to the convex combinations of intermediate points. + +\ccIncludeExampleCode{QP_solver/important_variables.cpp} + +Not surprisingly, only two points contribute to any convex +combination. \section{Working from MPS Files} diff --git a/QP_solver/examples/QP_solver/convex_hull_containment.cpp b/QP_solver/examples/QP_solver/convex_hull_containment.cpp index fc2b583c0fd..d169f78fe4d 100644 --- a/QP_solver/examples/QP_solver/convex_hull_containment.cpp +++ b/QP_solver/examples/QP_solver/convex_hull_containment.cpp @@ -3,11 +3,20 @@ #include #include #include -#include "is_in_convex_hull.h" +#include "solve_convex_hull_lp.h" typedef CGAL::Homogeneous_d Kernel_d; typedef Kernel_d::Point_d Point_d; +bool is_in_convex_hull (const Point_d& p, + std::vector::const_iterator begin, + std::vector::const_iterator end) +{ + CGAL::Quadratic_program_solution s = + solve_convex_hull_lp (p, begin, end, CGAL::MP_Float()); + return s.status() != CGAL::QP_INFEASIBLE; +} + int main() { std::vector points; @@ -19,7 +28,7 @@ int main() for (int j=0; j<=10; ++j) { // (i,j) is in the simplex iff i+j <= 10 bool contained = is_in_convex_hull - (Point_d (i, j), points.begin(), points.end(), CGAL::MP_Float()); + (Point_d (i, j), points.begin(), points.end()); assert (contained == (i+j<=10)); } diff --git a/QP_solver/examples/QP_solver/convex_hull_containment2.cpp b/QP_solver/examples/QP_solver/convex_hull_containment2.cpp index 234197d066c..fac6d24aa1b 100644 --- a/QP_solver/examples/QP_solver/convex_hull_containment2.cpp +++ b/QP_solver/examples/QP_solver/convex_hull_containment2.cpp @@ -3,11 +3,20 @@ #include #include #include -#include "is_in_convex_hull2.h" +#include "solve_convex_hull_lp2.h" typedef CGAL::Homogeneous_d Kernel_d; typedef Kernel_d::Point_d Point_d; +bool is_in_convex_hull (const Point_d& p, + std::vector::const_iterator begin, + std::vector::const_iterator end) +{ + CGAL::Quadratic_program_solution s = + solve_convex_hull_lp (p, begin, end, CGAL::MP_Float()); + return s.status() != CGAL::QP_INFEASIBLE; +} + int main() { std::vector points; @@ -19,7 +28,7 @@ int main() for (int j=0; j<=10; ++j) { // (i,j) is in the simplex iff i+j <= 10 bool contained = is_in_convex_hull - (Point_d (i, j), points.begin(), points.end(), CGAL::MP_Float()); + (Point_d (i, j), points.begin(), points.end()); assert (contained == (i+j<=10)); } diff --git a/QP_solver/examples/QP_solver/important_variables.cpp b/QP_solver/examples/QP_solver/important_variables.cpp new file mode 100644 index 00000000000..cb31ef121d8 --- /dev/null +++ b/QP_solver/examples/QP_solver/important_variables.cpp @@ -0,0 +1,37 @@ +// Example: find the points that contribute to a convex combination +#include +#include +#include +#include +#include "solve_convex_hull_lp2.h" + +typedef CGAL::Cartesian_d Kernel_d; +typedef Kernel_d::Point_d Point_d; +typedef CGAL::Quadratic_program_solution Solution; + +int main() +{ + std::vector points; + // convex hull: line spanned by {(0,0), (1,0),..., (9,0)} + for (int j=0; j<10; ++j) + points.push_back (Point_d (j, 0)); + + for (double f=0.5; f<10; ++f) { + Point_d p (f, 0.0); + Solution s = solve_convex_hull_lp + (p, points.begin(), points.end(), CGAL::MP_Float()); + std::cout << p; + if (s.status() == CGAL::QP_INFEASIBLE) + std::cout << " is not in the convex hull." << std::endl; + else { + std::cout << " is a convex combination of the points "; + Solution::Index_iterator it = s.basic_variable_indices_begin(); + Solution::Index_iterator end = s.basic_variable_indices_end(); + for (; it != end; ++it) + std::cout << points[*it] << " "; + std::cout << std::endl; + } + } + + return 0; +} diff --git a/QP_solver/examples/QP_solver/is_in_convex_hull.h b/QP_solver/examples/QP_solver/solve_convex_hull_lp.h similarity index 90% rename from QP_solver/examples/QP_solver/is_in_convex_hull.h rename to QP_solver/examples/QP_solver/solve_convex_hull_lp.h index cabe1217191..da6b2d793d0 100644 --- a/QP_solver/examples/QP_solver/is_in_convex_hull.h +++ b/QP_solver/examples/QP_solver/solve_convex_hull_lp.h @@ -36,8 +36,6 @@ solve_convex_hull_lp (const Point_d& p, typedef CGAL::Nonnegative_linear_program_from_iterators Program; - // the solution type - typedef CGAL::Quadratic_program_solution Solution; // ok, we are prepared now: construct program and solve it Program lp (end-begin, // number of variables @@ -46,7 +44,3 @@ solve_convex_hull_lp (const Point_d& p, R_it (CGAL::EQUAL), C_it (0)); return CGAL::solve_nonnegative_linear_program (lp, ET(0)); } - - // p is in the convex hull if and only if the program is feasible - return (s.status() != CGAL::QP_INFEASIBLE); -} diff --git a/QP_solver/examples/QP_solver/is_in_convex_hull2.h b/QP_solver/examples/QP_solver/solve_convex_hull_lp2.h similarity index 78% rename from QP_solver/examples/QP_solver/is_in_convex_hull2.h rename to QP_solver/examples/QP_solver/solve_convex_hull_lp2.h index 07f860e673c..20c2cf2b960 100644 --- a/QP_solver/examples/QP_solver/is_in_convex_hull2.h +++ b/QP_solver/examples/QP_solver/solve_convex_hull_lp2.h @@ -15,18 +15,16 @@ struct Homogeneous_begin { // function to test whether point is in the convex hull of other points; // the type ET is an exact type used for the computations template -bool is_in_convex_hull -(const Point_d& p, - RandomAccessIterator begin, RandomAccessIterator end, const ET& dummy) +CGAL::Quadratic_program_solution +solve_convex_hull_lp (const Point_d& p, + RandomAccessIterator begin, + RandomAccessIterator end, const ET& dummy) { // the right-hand side type typedef typename Point_d::Homogeneous_const_iterator B_it; - // the solution type - typedef CGAL::Quadratic_program_solution Solution; - // construct program and solve it - Solution s = CGAL::solve_nonnegative_linear_program + return CGAL::solve_nonnegative_linear_program (CGAL::make_nonnegative_linear_program_from_iterators (end-begin, // n p.dimension()+1, // m @@ -36,7 +34,5 @@ bool is_in_convex_hull CGAL::Const_oneset_iterator(CGAL::EQUAL), // ~ CGAL::Const_oneset_iterator ::value_type> (0)), ET(0)); // c - - // p is in the convex hull if and only if the program is feasible - return (s.status() != CGAL::QP_INFEASIBLE); } + diff --git a/QP_solver/include/CGAL/QP_solver.h b/QP_solver/include/CGAL/QP_solver.h index 79fbac23ddf..f81872a47f2 100644 --- a/QP_solver/include/CGAL/QP_solver.h +++ b/QP_solver/include/CGAL/QP_solver.h @@ -327,7 +327,7 @@ private: mutable Verbose_ostream vout2; // used for more diagnostic output mutable Verbose_ostream vout3; // used for full diagnostic output mutable Verbose_ostream vout4; // used for output of basis inverse - mutable Verbose_ostream vout5; // used for output of validity tests + mutable Verbose_ostream vout5; // used for output of validity tests // pricing strategy Pricing_strategy* strategyP; @@ -1277,9 +1277,9 @@ private: bool check_basis_inverse( Tag_false is_linear); // diagnostic output - void print_program ( ); - void print_basis ( ); - void print_solution( ); + void print_program ( ) const; + void print_basis ( ) const; + void print_solution( ) const; void print_ratio_1_original(int k, const ET& x_k, const ET& q_k); void print_ratio_1_slack(int k, const ET& x_k, const ET& q_k); diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h index 8e1e16c209a..613271f752b 100644 --- a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h +++ b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h @@ -2972,7 +2972,7 @@ set_verbosity( int verbose, std::ostream& stream) template < typename Q, typename ET, typename Tags > void QP_solver:: -print_program( ) +print_program( ) const { int row, i; @@ -3050,7 +3050,7 @@ print_program( ) template < typename Q, typename ET, typename Tags > void QP_solver:: -print_basis( ) +print_basis( ) const { char label; vout1 << " basis: "; @@ -3069,7 +3069,8 @@ print_basis( ) if ( has_ineq) { vout2.out() << std::endl << "basic constraints: "; - for (Index_iterator i_it = C.begin(); i_it != C.end(); ++i_it) { + for (Index_const_iterator i_it = + C.begin(); i_it != C.end(); ++i_it) { label = (qp_r[*i_it] == CGAL::EQUAL) ? 'e' : 'i'; vout2.out() << *i_it << ":" << label << " "; } @@ -3103,7 +3104,7 @@ print_basis( ) template < typename Q, typename ET, typename Tags > void QP_solver:: -print_solution( ) +print_solution( ) const { if ( vout3.verbose()) { vout3.out() << std::endl