diff --git a/.gitattributes b/.gitattributes index 540003b0f20..97821d59308 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1681,6 +1681,10 @@ Principal_component_analysis/demo/Principal_component_analysis/windows/3d/res/id Principal_component_analysis/demo/Principal_component_analysis/windows/3d/res/idb8.bmp -text svneol=unset#image/bmp Principal_component_analysis/doc_tex/Principal_component_analysis/fit.eps -text svneol=unset#application/postscript Principal_component_analysis/doc_tex/Principal_component_analysis/fit.png -text svneol=unset#image/png +QP_solver/doc_tex/QP_solver/closest_point.eps -text +QP_solver/doc_tex/QP_solver/closest_point.fig -text +QP_solver/doc_tex/QP_solver/closest_point.gif -text +QP_solver/doc_tex/QP_solver/closest_point.pdf -text QP_solver/doc_tex/QP_solver/first_lp.eps -text QP_solver/doc_tex/QP_solver/first_lp.fig -text QP_solver/doc_tex/QP_solver/first_lp.gif -text @@ -1713,6 +1717,7 @@ 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/rational_qp_solver.cin -text QP_solver/examples/QP_solver/rational_qp_solver.data -text +QP_solver/examples/QP_solver/shortest_vector.cpp -text QP_solver/examples/QP_solver/solve_convex_hull_containment_lp.h -text QP_solver/examples/QP_solver/solve_convex_hull_containment_lp2.h -text QP_solver/include/CGAL/QP_functions.h -text diff --git a/QP_solver/doc_tex/QP_solver/closest_point.eps b/QP_solver/doc_tex/QP_solver/closest_point.eps new file mode 100644 index 00000000000..48cff8c50e3 --- /dev/null +++ b/QP_solver/doc_tex/QP_solver/closest_point.eps @@ -0,0 +1,186 @@ +%!PS-Adobe-2.0 EPSF-2.0 +%%Title: closest_point.fig +%%Creator: fig2dev Version 3.2 Patchlevel 4 +%%CreationDate: Fri Oct 27 16:03:35 2006 +%%For: gaertner@tamtam.inf.ethz.ch (Bernd Gaertner) +%%BoundingBox: 0 0 326 218 +%%Magnification: 1.0000 +%%EndComments +/$F2psDict 200 dict def +$F2psDict begin +$F2psDict /mtrx matrix put +/col-1 {0 setgray} bind def +/col0 {0.000 0.000 0.000 srgb} bind def +/col1 {0.000 0.000 1.000 srgb} bind def +/col2 {0.000 1.000 0.000 srgb} bind def +/col3 {0.000 1.000 1.000 srgb} bind def +/col4 {1.000 0.000 0.000 srgb} bind def +/col5 {1.000 0.000 1.000 srgb} bind def +/col6 {1.000 1.000 0.000 srgb} bind def +/col7 {1.000 1.000 1.000 srgb} bind def +/col8 {0.000 0.000 0.560 srgb} bind def +/col9 {0.000 0.000 0.690 srgb} bind def +/col10 {0.000 0.000 0.820 srgb} bind def +/col11 {0.530 0.810 1.000 srgb} bind def +/col12 {0.000 0.560 0.000 srgb} bind def +/col13 {0.000 0.690 0.000 srgb} bind def +/col14 {0.000 0.820 0.000 srgb} bind def +/col15 {0.000 0.560 0.560 srgb} bind def +/col16 {0.000 0.690 0.690 srgb} bind def +/col17 {0.000 0.820 0.820 srgb} bind def +/col18 {0.560 0.000 0.000 srgb} bind def +/col19 {0.690 0.000 0.000 srgb} bind def +/col20 {0.820 0.000 0.000 srgb} bind def +/col21 {0.560 0.000 0.560 srgb} bind def +/col22 {0.690 0.000 0.690 srgb} bind def +/col23 {0.820 0.000 0.820 srgb} bind def +/col24 {0.500 0.190 0.000 srgb} bind def +/col25 {0.630 0.250 0.000 srgb} bind def +/col26 {0.750 0.380 0.000 srgb} bind def +/col27 {1.000 0.500 0.500 srgb} bind def +/col28 {1.000 0.630 0.630 srgb} bind def +/col29 {1.000 0.750 0.750 srgb} bind def +/col30 {1.000 0.880 0.880 srgb} bind def +/col31 {1.000 0.840 0.000 srgb} bind def + +end +save +newpath 0 218 moveto 0 0 lineto 326 0 lineto 326 218 lineto closepath clip newpath +-107.3 378.7 translate +1 -1 scale + +/cp {closepath} bind def +/ef {eofill} bind def +/gr {grestore} bind def +/gs {gsave} bind def +/sa {save} bind def +/rs {restore} bind def +/l {lineto} bind def +/m {moveto} bind def +/rm {rmoveto} bind def +/n {newpath} bind def +/s {stroke} bind def +/sh {show} bind def +/slc {setlinecap} bind def +/slj {setlinejoin} bind def +/slw {setlinewidth} bind def +/srgb {setrgbcolor} bind def +/rot {rotate} bind def +/sc {scale} bind def +/sd {setdash} bind def +/ff {findfont} bind def +/sf {setfont} bind def +/scf {scalefont} bind def +/sw {stringwidth} bind def +/tr {translate} bind def +/tnt {dup dup currentrgbcolor + 4 -2 roll dup 1 exch sub 3 -1 roll mul add + 4 -2 roll dup 1 exch sub 3 -1 roll mul add + 4 -2 roll dup 1 exch sub 3 -1 roll mul add srgb} + bind def +/shd {dup dup currentrgbcolor 4 -2 roll mul 4 -2 roll mul + 4 -2 roll mul srgb} bind def + /DrawEllipse { + /endangle exch def + /startangle exch def + /yrad exch def + /xrad exch def + /y exch def + /x exch def + /savematrix mtrx currentmatrix def + x y tr xrad yrad sc 0 0 1 startangle endangle arc + closepath + savematrix setmatrix + } def + +/$F2psBegin {$F2psDict begin /$F2psEnteredState save def} def +/$F2psEnd {$F2psEnteredState restore end} def + +$F2psBegin +10 setmiterlimit +0 slj 0 slc + 0.06000 0.06000 sc +% +% Fig objects follow +% +% +% here starts figure with depth 50 +% Polyline +7.500 slw +n 2400 2700 m + 2625 2700 l gs col0 s gr +% Polyline +n 2400 2775 m + 2625 2775 l gs col0 s gr +% Polyline +n 2400 2850 m + 2625 2850 l gs col0 s gr +% Polyline +n 7200 5475 m + 7200 5700 l gs col0 s gr +% Polyline +n 7125 5475 m + 7125 5700 l gs col0 s gr +% Polyline +n 7050 5475 m + 7050 5700 l gs col0 s gr +% Polyline +n 5700 2850 m + 5775 3000 l gs col0 s gr +% Polyline +n 5850 2775 m + 5925 2925 l gs col0 s gr +% Polyline +n 6000 2700 m + 6075 2850 l gs col0 s gr +% Polyline +n 1800 5700 m + 7200 5700 l gs col0 s gr +% Polyline +n 2400 2700 m + 2400 6300 l gs col0 s gr +% Polyline +n 1800 4800 m + 6000 2700 l gs col0 s gr +% Polyline +n 7200 6300 m + 3600 2700 l gs col0 s gr +% Polyline +n 3600 2700 m + 3450 2850 l gs col0 s gr +% Polyline +n 3675 2775 m + 3525 2925 l gs col0 s gr +% Polyline +n 3750 2850 m + 3600 3000 l gs col0 s gr +% Polyline +n 2400 5700 m 2400 4500 l 4395 3495 l 6600 5700 l + cp gs col12 0.25 tnt ef gr gs col0 s gr +% Polyline +15.000 slw + [90] 0 sd +n 6600 3900 m + 5700 4800 l gs col0 s gr [] 0 sd +/Times-Italic ff 300.00 scf sf +5925 4875 m +gs 1 -1 sc (q) col0 sh gr +/Times-Italic ff 300.00 scf sf +6825 4125 m +gs 1 -1 sc (p) col0 sh gr +% here ends figure; +% +% here starts figure with depth 47 +% Ellipse +7.500 slw + [60] 0 sd +n 6600 3900 75 75 0 360 DrawEllipse gs 0.00 setgray ef gr gs col0 s gr + [] 0 sd +% Ellipse + [60] 0 sd +n 5700 4800 75 75 0 360 DrawEllipse gs 0.00 setgray ef gr gs col0 s gr + [] 0 sd +% here ends figure; +$F2psEnd +rs +showpage diff --git a/QP_solver/doc_tex/QP_solver/closest_point.fig b/QP_solver/doc_tex/QP_solver/closest_point.fig new file mode 100644 index 00000000000..fa4ec1dfb55 --- /dev/null +++ b/QP_solver/doc_tex/QP_solver/closest_point.fig @@ -0,0 +1,55 @@ +#FIG 3.2 +Portrait +Center +Inches +Letter +100.00 +Single +-2 +1200 2 +6 2400 2700 2625 2850 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 2400 2700 2625 2700 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 2400 2775 2625 2775 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 2400 2850 2625 2850 +-6 +6 7050 5475 7200 5700 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 7200 5475 7200 5700 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 7125 5475 7125 5700 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 7050 5475 7050 5700 +-6 +6 5700 2700 6075 3000 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 5700 2850 5775 3000 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 5850 2775 5925 2925 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 6000 2700 6075 2850 +-6 +1 3 1 1 0 0 47 -1 20 4.000 1 0.0000 6600 3900 75 75 6600 3900 6675 3900 +1 3 1 1 0 0 47 -1 20 4.000 1 0.0000 5700 4800 75 75 5700 4800 5775 4800 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 1800 5700 7200 5700 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 2400 2700 2400 6300 +2 1 0 1 0 7 50 -1 -1 4.000 0 0 -1 0 0 2 + 1800 4800 6000 2700 +2 1 0 1 0 0 50 -1 -1 4.000 0 0 -1 0 0 2 + 7200 6300 3600 2700 +2 1 0 1 0 0 50 -1 -1 4.000 0 0 -1 0 0 2 + 3600 2700 3450 2850 +2 1 0 1 0 0 50 -1 -1 4.000 0 0 -1 0 0 2 + 3675 2775 3525 2925 +2 1 0 1 0 0 50 -1 -1 4.000 0 0 -1 0 0 2 + 3750 2850 3600 3000 +2 3 0 1 0 12 50 -1 25 4.000 0 0 -1 0 0 5 + 2400 5700 2400 4500 4395 3495 6600 5700 2400 5700 +2 1 1 2 0 7 50 -1 -1 6.000 0 0 -1 0 0 2 + 6600 3900 5700 4800 +4 0 0 50 -1 1 20 0.0000 4 120 90 5925 4875 q\001 +4 0 0 50 -1 1 20 0.0000 4 120 90 6825 4125 p\001 diff --git a/QP_solver/doc_tex/QP_solver/closest_point.gif b/QP_solver/doc_tex/QP_solver/closest_point.gif new file mode 100644 index 00000000000..c01f26a12c0 Binary files /dev/null and b/QP_solver/doc_tex/QP_solver/closest_point.gif differ diff --git a/QP_solver/doc_tex/QP_solver/closest_point.pdf b/QP_solver/doc_tex/QP_solver/closest_point.pdf new file mode 100644 index 00000000000..fb6f7a07088 Binary files /dev/null and b/QP_solver/doc_tex/QP_solver/closest_point.pdf differ diff --git a/QP_solver/doc_tex/QP_solver/main.tex b/QP_solver/doc_tex/QP_solver/main.tex index 84d55791962..f2484e82fdd 100644 --- a/QP_solver/doc_tex/QP_solver/main.tex +++ b/QP_solver/doc_tex/QP_solver/main.tex @@ -53,13 +53,21 @@ the smallest objective function value among all feasible solutions. Such a solution may not exist, see Sections \ref{sec:QP-infeasible} and \ref{sec:QP-unbounded} below. -\subsection{Linear and nonnegative programs.} +\subsection{Linear, nonnegative, and free programs.} If $D=0$, the quadratic program is in fact a \emph{linear program}, and in the case that $l$ is the zero vector and all entries of $u$ are $\infty$, the program is said to be \emph{nonnegative}. The package offers dedicated solution methods for these special cases, see Sections \ref{sec:QP-lp} and \ref{sec:QP-nonnegative} below. +Often, there are no bounds at all for the variables, i.e.\ all entries +of $l$ are $-\infty$, and all entries of $u$ are $\infty$. Such a +program is called \emph{free}. There is no dedicated solution method +for this case (a free quadratic or linear program is treated like a +general quadratic or linear program), but the package provides models +for free programs that save you from specifying the infinite bounds (see +Section \ref{sec:QP-important_constraints} for an example). + \subsection{A first example} Let's consider the following quadratic program in two variables: \[ @@ -379,6 +387,32 @@ contribute to the convex combinations of intermediate points. Not surprisingly, only two points contribute to any convex combination. +\section{The Important Constraints}\label{sec:QP-important_constraints} +If we have a linear or quadratic program with many inequality constraints, +the ``important'' constraints are typically the ones that are satisfied with +equality at the optimal solution. + +To demonstrate this, let's consider the problem of finding the point $q$ +in the intersection of halfspaces that is closest to a given point $p$, +see Figure \ref{fig:QP-closest_point}. + +\begin{figure}[htbp] +\begin{ccTexOnly} +\begin{center} +\includegraphics{QP_solver/closest_point} +\end{center} +\end{ccTexOnly} +\caption{The closest point in the intersection of halfspaces +\label{fig:QP-closest_point}} + +\begin{ccHtmlOnly} +
+The closest point in the intersection of halfspaces +
+\end{ccHtmlOnly} +\end{figure} + + \section{Working from MPS Files} \section{Infeasibility}\label{sec:QP-infeasible} @@ -386,7 +420,7 @@ combination. \section{Unboundedness}\label{sec:QP-unbounded} -\section{The Important Constraints} + \section{Solution Certificates} diff --git a/QP_solver/examples/QP_solver/shortest_vector.cpp b/QP_solver/examples/QP_solver/shortest_vector.cpp new file mode 100644 index 00000000000..250f8779811 --- /dev/null +++ b/QP_solver/examples/QP_solver/shortest_vector.cpp @@ -0,0 +1,90 @@ +// Example: computes point of minimum norm in the intersection of halfspaces +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +// unary function that maps integer j to iterator constructed from j +template +struct Jth { + typedef Iterator result_type; + result_type operator() (int j) const { return Iterator (j); } +}; + +// unary function that maps a hyperplane h to its negative j-th coordinate, +// for some fixed j +template +class Fixed_coordinate { +public: + typedef typename CGAL::Kernel_traits::Kernel::RT result_type; + Fixed_coordinate (int index) : j (index) {} + result_type operator() (const Hyperplane& h) const {return -h[j];} +private: + const int j; +}; + +// unary function that maps an integer i to 1 if i=j and 0 otherwise, +// for some fixed j +template +class Unit_vector { +public: + typedef RT result_type; + Unit_vector (int index) : j (index) {} + result_type operator() (int i) { return i==j ? 1 : 0; } +private: + const int j; +}; + +// function to solve the QP that computes the minimum-norm point in the +// intersection of positive halfspaces induced by a set of oriented +// hyperplanes +template +CGAL::Quadratic_program_solution +solve_minimum_norm_point_qp +(HyperplaneIterator begin, HyperplaneIterator end, const ET& dummy) +{ + // hyperplane type + typedef typename HyperplaneIterator::value_type H; + // input number type + typedef typename CGAL::Kernel_traits::Kernel::RT RT; + // iterator type for j-th column of A (j=0,...,d-1) and for b (j=d); + // transforms iterator for hyperplane to iterator for j-th coordinates + typedef CGAL::Join_input_iterator_1 + > Fixed_coordinate_iterator; + // iterator type for A; transforms iterator for indices to iterator + // for corresponding columns of A + typedef CGAL::Join_input_iterator_1 + , Jth > A_it; + // iterator type for b; + typedef Fixed_coordinate_iterator B_it; + // iterator type for relations (all are "<=") + typedef CGAL::Const_oneset_iterator R_it; + // iterator type for c (all entries are 0) + typedef CGAL::Const_oneset_iterator C_it; + // iterator type for i-th row of D; transforms iterator for indices + // to iterator for the i-th unit vector + typedef CGAL::Join_input_iterator_1 + , Unit_vector > Unit_vector_iterator; + // iterator type for D; transforms iterator for indices to iterator + // for corresponding rows of D + typedef CGAL::Join_input_iterator_1 + , Jth > D_it; + // determine dimension from first hyperplane + int d = (*begin).dimension(); + // construct program and solve itx + return CGAL::solve_quadratic_program + (CGAL::make_free_quadratic_program_from_iterators + (d, end-begin, A_it(0), B_it(d), R_it(CGAL::SMALLER), D_it(), C_it (0))); +} + +int main() { + // + return 0; + +} diff --git a/QP_solver/include/CGAL/QP_models.h b/QP_solver/include/CGAL/QP_models.h index f1f3d487310..fc27e6f679d 100644 --- a/QP_solver/include/CGAL/QP_models.h +++ b/QP_solver/include/CGAL/QP_models.h @@ -349,6 +349,7 @@ public: d, c, c0) {} }; + // corresponding global function // make_nonnegative_quadratic_program_from_iterators // ------------------------------------------------- @@ -400,6 +401,98 @@ public: {} }; +// Free_quadratic_program_from_iterators +// ------------------------------------- +template < + typename A_it, // for constraint matrix A (columnwise) + typename B_it, // for right-hand side b + typename R_it, // for relations (value type Comparison) + typename D_it, // for quadratic matrix D (rowwise) + typename C_it > // for objective function c +class Free_quadratic_program_from_iterators : + public Quadratic_program_from_iterators ::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + D_it, C_it> +{ +private: + typedef Quadratic_program_from_iterators ::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + D_it, C_it> Base; + typedef typename QP_model_default_iterators::It_1d Const_FLU_iterator; + typedef typename QP_model_default_iterators::It_1d Const_LU_iterator; +public: + QP_MODEL_ITERATOR_TYPES; + Free_quadratic_program_from_iterators ( + int n, int m, // number of variables / constraints + const A_iterator& a, + const B_iterator& b, + const R_iterator& r, + const D_iterator& d, + const C_iterator& c, + const C_entry& c0 = C_entry(0) + ) + : Base (n, m, a, b, r, + Const_FLU_iterator(false), Const_LU_iterator(C_entry(0)), + Const_FLU_iterator(false), Const_LU_iterator(C_entry(0)), + d, c, c0) + {} +}; +// corresponding global function +// make_free_quadratic_program_from_iterators +// ------------------------------------------ +template < + typename A_it, // for constraint matrix A (columnwise) + typename B_it, // for right-hand side b + typename R_it, // for relations (value type Comparison) + typename D_it, // for quadratic matrix D (rowwise) + typename C_it > // for objective function c +Free_quadratic_program_from_iterators + +make_free_quadratic_program_from_iterators ( + int n, int m, + const A_it& a, + const B_it& b, + const R_it& r, + const D_it& d, + const C_it& c, + typename std::iterator_traits::value_type c0 = + typename std::iterator_traits::value_type(0)) +{ + return Free_quadratic_program_from_iterators + + (n, m, a, b, r, d, c, c0); +} + +// Free_Quadratic_program_from_pointers +// ------------------------------------------- +template +class Free_quadratic_program_from_pointers : + public Free_quadratic_program_from_iterators + +{ +private: + typedef Free_quadratic_program_from_iterators + Base; +public: + QP_MODEL_ITERATOR_TYPES; + Free_quadratic_program_from_pointers ( + int n, int m, // number of variables / constraints + const A_iterator& a, + const B_iterator& b, + const R_iterator& r, + const D_iterator& d, + const C_iterator& c, + const C_entry& c0 = C_entry(0) + ) + : Base (n, m, a, b, r, d, c, c0) + {} +}; // Nonnegative_linear_program_from_iterators @@ -492,6 +585,96 @@ public: {} }; +// Free_linear_program_from_iterators +// ---------------------------------- +template < + typename A_it, // for constraint matrix A (columnwise) + typename B_it, // for right-hand side b + typename R_it, // for relations (value type Comparison) + typename C_it > // for objective function c +class Free_linear_program_from_iterators : + public Quadratic_program_from_iterators ::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_2d, C_it> +{ +private: + typedef Quadratic_program_from_iterators ::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_1d, + typename QP_model_default_iterators::It_2d, C_it> Base; + typedef typename QP_model_default_iterators::It_1d Const_FLU_iterator; + typedef typename QP_model_default_iterators::It_1d Const_LU_iterator; + typedef typename QP_model_default_iterators::It_2d Const_D_iterator; +public: + QP_MODEL_ITERATOR_TYPES; + Free_linear_program_from_iterators ( + int n, int m, // number of variables / constraints + const A_iterator& a, + const B_iterator& b, + const R_iterator& r, + const C_iterator& c, + const C_entry& c0 = C_entry(0) + ) + : Base (n, m, a, b, r, + Const_FLU_iterator(false), Const_LU_iterator(C_entry(0)), + Const_FLU_iterator(false), Const_LU_iterator(C_entry(0)), + Const_D_iterator(C_entry(0)), c, c0) + {} +}; + +// corresponding global function +// make_free_linear_program_from_iterators +// --------------------------------------- +template < + typename A_it, // for constraint matrix A (columnwise) + typename B_it, // for right-hand side b + typename R_it, // for relations (value type Comparison) + typename C_it > // for objective function c +Nonnegative_linear_program_from_iterators + +make_free_linear_program_from_iterators ( + int n, int m, + const A_it& a, + const B_it& b, + const R_it& r, + const C_it& c, + typename std::iterator_traits::value_type c0 = + typename std::iterator_traits::value_type(0)) +{ + return Free_linear_program_from_iterators + + (n, m, a, b, r, c, c0); +} + +// Free_linear_program_from_pointers +// --------------------------------- +template +class Free_linear_program_from_pointers : + public Free_linear_program_from_iterators + +{ +private: + typedef Free_linear_program_from_iterators + Base; +public: + QP_MODEL_ITERATOR_TYPES; + Free_linear_program_from_pointers ( + int n, int m, // number of variables / constraints + const A_iterator& a, + const B_iterator& b, + const R_iterator& r, + const C_iterator& c, + const C_entry& c0 = C_entry(0) + ) + : Base (n, m, a, b, r, c, c0) + {} +}; + // Quadratic_program_from_mps // -------------------------- namespace QP_from_mps_detail {