- Free models added;

- figures/program for closest point in intersection of halfspaces
This commit is contained in:
Bernd Gärtner 2006-10-27 14:07:45 +00:00
parent 3a8a8f00d2
commit aeb1213f4c
8 changed files with 555 additions and 2 deletions

5
.gitattributes vendored
View File

@ -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

View File

@ -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

View File

@ -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

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.7 KiB

Binary file not shown.

View File

@ -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}
<CENTER>
<IMG BORDER=0 SRC="./closest_point.gif" ALIGN=center ALT="The closest point in the intersection of halfspaces">
</CENTER>
\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}

View File

@ -0,0 +1,90 @@
// Example: computes point of minimum norm in the intersection of halfspaces
#include <cassert>
#include <vector>
#include <CGAL/Cartesian_d.h>
#include <CGAL/MP_Float.h>
#include <boost/iterator/counting_iterator.hpp>
#include <CGAL/iterator.h>
#include <CGAL/Kernel_traits.h>
#include <CGAL/QP_models.h>
#include <CGAL/QP_functions.h>
// unary function that maps integer j to iterator constructed from j
template <class Iterator>
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 Hyperplane>
class Fixed_coordinate {
public:
typedef typename CGAL::Kernel_traits<Hyperplane>::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 RT>
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 <class HyperplaneIterator, class ET>
CGAL::Quadratic_program_solution<ET>
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<H>::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
<HyperplaneIterator, Fixed_coordinate<H> > Fixed_coordinate_iterator;
// iterator type for A; transforms iterator for indices to iterator
// for corresponding columns of A
typedef CGAL::Join_input_iterator_1
<boost::counting_iterator<int>, Jth<Fixed_coordinate_iterator> > A_it;
// iterator type for b;
typedef Fixed_coordinate_iterator B_it;
// iterator type for relations (all are "<=")
typedef CGAL::Const_oneset_iterator<CGAL::Comparison_result> R_it;
// iterator type for c (all entries are 0)
typedef CGAL::Const_oneset_iterator<RT> 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
<boost::counting_iterator<int>, Unit_vector<RT> > Unit_vector_iterator;
// iterator type for D; transforms iterator for indices to iterator
// for corresponding rows of D
typedef CGAL::Join_input_iterator_1
<boost::counting_iterator<int>, Jth<Unit_vector_iterator> > 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;
}

View File

@ -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 <A_it, B_it, R_it,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
D_it, C_it>
{
private:
typedef Quadratic_program_from_iterators <A_it, B_it, R_it,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
D_it, C_it> Base;
typedef typename QP_model_default_iterators<bool*>::It_1d Const_FLU_iterator;
typedef typename QP_model_default_iterators<B_it>::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
<A_it, B_it, R_it, D_it, C_it>
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<C_it>::value_type c0 =
typename std::iterator_traits<C_it>::value_type(0))
{
return Free_quadratic_program_from_iterators
<A_it, B_it, R_it, D_it, C_it>
(n, m, a, b, r, d, c, c0);
}
// Free_Quadratic_program_from_pointers
// -------------------------------------------
template <typename NT>
class Free_quadratic_program_from_pointers :
public Free_quadratic_program_from_iterators
<NT**, NT*, CGAL::Comparison_result*, NT**, NT*>
{
private:
typedef Free_quadratic_program_from_iterators
<NT**, NT*, CGAL::Comparison_result*, NT**, NT*> 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 <A_it, B_it, R_it,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<B_it>::It_2d, C_it>
{
private:
typedef Quadratic_program_from_iterators <A_it, B_it, R_it,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<bool*>::It_1d,
typename QP_model_default_iterators<B_it>::It_1d,
typename QP_model_default_iterators<B_it>::It_2d, C_it> Base;
typedef typename QP_model_default_iterators<bool*>::It_1d Const_FLU_iterator;
typedef typename QP_model_default_iterators<B_it>::It_1d Const_LU_iterator;
typedef typename QP_model_default_iterators<B_it>::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
<A_it, B_it, R_it, C_it>
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<C_it>::value_type c0 =
typename std::iterator_traits<C_it>::value_type(0))
{
return Free_linear_program_from_iterators
<A_it, B_it, R_it, C_it>
(n, m, a, b, r, c, c0);
}
// Free_linear_program_from_pointers
// ---------------------------------
template <typename NT>
class Free_linear_program_from_pointers :
public Free_linear_program_from_iterators
<NT**, NT*, CGAL::Comparison_result*, NT*>
{
private:
typedef Free_linear_program_from_iterators
<NT**, NT*, CGAL::Comparison_result*, NT*> 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 {