songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5741 from danston/Solvers-add_osqp-danston 

[Small Feature] OSQP Support in Solver Interface
This commit is contained in:
Laurent Rineau 2021-07-27 15:55:23 +02:00
parent 249d0f3af0 5e17827c64
commit b9743fffa3
30 changed files with 796 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
Documentation/doc/Documentation/Third_party.txt
View File

@ -297,11 +297,8 @@ If you experience such an issue, we recommand to compile \ceres without `glog` s
\glpk (GNU Linear Programming Kit) is a library for solving linear programming (LP), mixed integer programming (MIP), and other related problems.
In \cgal, \glpk provides an optional linear integer program solver
in the \ref PkgPolygonalSurfaceReconstruction package. In order to use
\glpk in \cgal programs, the executables should be linked with the
CMake imported target `CGAL::GLPK_support` provided in
`CGAL_GLPK_support.cmake`.
In \cgal, \glpk provides an optional linear integer program solver in the \ref PkgPolygonalSurfaceReconstruction package.
In order to use \glpk in \cgal programs, the executables should be linked with the CMake imported target `CGAL::GLPK_support` provided in `CGAL_GLPK_support.cmake`.
The \glpk web site is <A HREF="https://www.gnu.org/software/glpk/">`https://www.gnu.org/software/glpk/`</A>.
@ -309,12 +306,18 @@ The \glpk web site is <A HREF="https://www.gnu.org/software/glpk/">`https://www.
\scip (Solving Constraint Integer Programs) is currently one of the fastest open source solvers for mixed integer programming (MIP) and mixed integer nonlinear programming (MINLP).
In \cgal, \scip provides an optional linear integer program solver
in the \ref PkgPolygonalSurfaceReconstruction package. In order to use
\scip in \cgal programs, the executables should be linked with the
CMake imported target `CGAL::SCIP_support` provided in
`CGAL_SCIP_support.cmake`.
In \cgal, \scip provides an optional linear integer program solver in the \ref PkgPolygonalSurfaceReconstruction package.
In order to use \scip in \cgal programs, the executables should be linked with the CMake imported target `CGAL::SCIP_support` provided in `CGAL_SCIP_support.cmake`.
The \scip web site is <A HREF="http://scip.zib.de/">`http://scip.zib.de/`</A>.
\subsection thirdpartyOSQP OSQP
\osqp (Operator Splitting Quadratic Program) is currently one of the fastest open source solvers for convex Quadratic Programs (QP).
In \cgal, \osqp provides an optional solver for the QP problems often arising in various computational geometry algorithms.
In order to use \osqp in \cgal programs, the executables should be linked with the CMake imported target `CGAL::OSQP_support` provided in `CGAL_OSQP_support.cmake`.
The \osqp web site is <A HREF="https://osqp.org">`https://osqp.org`</A>.
*/

1
Documentation/doc/resources/1.8.13/BaseDoxyfile.in
View File

@ -267,6 +267,7 @@ ALIASES = "cgal=%CGAL" \
"ceres=Ceres" \
"glpk=GLPK" \
"scip=SCIP" \
"osqp=OSQP" \
"rs=RS" \
"rs3=RS3" \
"unix=Unix" \

1
Documentation/doc/resources/1.8.14/BaseDoxyfile.in
View File

@ -268,6 +268,7 @@ ALIASES = "cgal=%CGAL" \
"ceres=Ceres" \
"glpk=GLPK" \
"scip=SCIP" \
"osqp=OSQP" \
"rs=RS" \
"rs3=RS3" \
"unix=Unix" \

1
Documentation/doc/resources/1.8.20/BaseDoxyfile.in
View File

@ -290,6 +290,7 @@ ALIASES = "cgal=%CGAL" \
"ceres=Ceres" \
"glpk=GLPK" \
"scip=SCIP" \
"osqp=OSQP" \
"rs=RS" \
"rs3=RS3" \
"unix=Unix" \

1
Documentation/doc/resources/1.8.4/BaseDoxyfile.in
View File

@ -236,6 +236,7 @@ ALIASES += "zlib=zlib"
ALIASES += "ceres=Ceres"
ALIASES += "glpk=GLPK"
ALIASES += "scip=SCIP"
ALIASES += "osqp=OSQP"
ALIASES += "rs=RS"
ALIASES += "rs3=RS3"
ALIASES += "unix=Unix"

5
Installation/CHANGES.md
View File

@ -9,6 +9,11 @@ Release date: December 2021
- Added the function `CGAL::Polygon_mesh_processing::match_faces()`, which, given two polygon meshes, identifies their common faces as well as as faces present in only either of them.
### [CGAL and Solvers](https://doc.cgal.org/5.4/Manual/packages.html#PkgSolverInterface)
- Added the [OSQP solver](https://osqp.org/) support. This solver enables to efficiently compute the convex Quadratic Programming (QP) problems arising in the context of several packages.
[Release 5.3](https://github.com/CGAL/cgal/releases/tag/v5.3)
-----------

7
Installation/cmake/modules/CGAL_OSQP_support.cmake Normal file
View File

@ -0,0 +1,7 @@
if(OSQP_FOUND AND NOT TARGET CGAL::OSQP_support)
add_library(CGAL::OSQP_support INTERFACE IMPORTED)
set_target_properties(CGAL::OSQP_support PROPERTIES
INTERFACE_COMPILE_DEFINITIONS "CGAL_USE_OSQP"
INTERFACE_INCLUDE_DIRECTORIES "${OSQP_INCLUDE_DIR}"
INTERFACE_LINK_LIBRARIES "${OSQP_LIBRARIES}")
endif()

23
Installation/cmake/modules/FindOSQP.cmake Normal file
View File

@ -0,0 +1,23 @@
# This file sets up OSQP for CMake. Once done this will define
# OSQP_FOUND - system has OSQP lib
# OSQP_INCLUDE_DIR - the OSQP include directory
# OSQP_LIBRARIES - link these to use OSQP
if(NOT OSQP_FOUND)
find_path(OSQP_INCLUDE_DIR
NAMES osqp.h
PATHS /usr/local/include/osqp/
ENV OSQP_INC_DIR)
find_library(OSQP_LIBRARIES
NAMES libosqp osqp
PATHS ENV LD_LIBRARY_PATH
ENV LIBRARY_PATH
/usr/local/lib
${OSQP_INCLUDE_DIR}/../lib
ENV OSQP_LIB_DIR)
if(OSQP_LIBRARIES AND OSQP_INCLUDE_DIR)
set(OSQP_FOUND TRUE)
endif()
endif()

91
Solver_interface/doc/Solver_interface/Concepts/LinearProgramTraits.h Normal file
View File

@ -0,0 +1,91 @@
/*!
\ingroup PkgSolverInterfaceConcepts
\cgalConcept
A concept that describes the set of methods used to define and solve a
linear programming (`lp`) problem of the general form:
<center>
\f{eqnarray*}{
& \mbox{minimize} & \mathbf{q}^{T}\mathbf{x} + r \\
& \mbox{subject to} & \mathbf{l} \leq A\mathbf{x} \leq \mathbf{u}
\f}
</center>
in \f$ n \f$ real variables \f$ \mathbf{x} = (x_0, \ldots, x_{n-1}) \f$ and \f$ m \f$ constraints.
Here,
<UL>
<LI>\f$ \mathbf{q} \f$ is an \f$ n \f$-dimensional vector (the linear objective function),
<LI>\f$ r \f$ is a constant,
<LI>\f$ A \f$ is an \f$ m\times n\f$ matrix (the constraint matrix),
<LI>\f$ \mathbf{l} \f$ is an \f$ m \f$-dimensional vector of lower constraint bounds,
where \f$ l_i \in \mathbb{R} \cup \{-\infty\} \f$ for all \f$ i \f$,
<LI>\f$ \mathbf{u} \f$ is an \f$ m \f$-dimensional vector of upper constraint bounds,
where \f$ u_i \in \mathbb{R} \cup \{+\infty\} \f$ for all \f$ i \f$.
</UL>
*/
class LinearProgramTraits {
public:
/// \name Memory
/// @{
/*!
Allocates memory for `n` variables and `m` constraints in `lp`.
*/
void resize(const std::size_t n, const std::size_t m) { }
/// @}
/// \name Initialization
/// @{
/*!
Sets the entry `qi` of `lp` to `value`.
*/
void set_q(const std::size_t i, const FT value) { }
/*!
Sets the entry `r` of `lp` to `value`.
*/
void set_r(const FT value) { }
/*!
Sets the entry `Aij` in the row `i` and column `j` of
the constraint matrix `A` of `lp` to `value`.
*/
void set_A(const std::size_t i, const std::size_t j, const FT value) { }
/*!
Sets the entry `li` of `lp` to `value`.
*/
void set_l(const std::size_t i, const FT value) { }
/*!
Sets the entry `ui` of `lp` to `value`.
*/
void set_u(const std::size_t i, const FT value) { }
/// @}
/// \name Solution
/// @{
/*!
\brief solves the linear program.
Number of values in `solution` equals to the number `n` of values in the vector `x`.
\tparam OutIterator
a model of `OutputIterator` that accepts values of type `FieldNumberType`
\param solution
an output iterator with the solution
\returns a status of the computation `success == true`
*/
template<typename OutIterator>
bool solve(OutIterator solution) { }
/// @}
};

100
Solver_interface/doc/Solver_interface/Concepts/QuadraticProgramTraits.h Normal file
View File

@ -0,0 +1,100 @@
/*!
\ingroup PkgSolverInterfaceConcepts
\cgalConcept
A concept that describes the set of methods used to define and solve a
quadratic programming (`qp`) problem of the general form:
<center>
\f{eqnarray*}{
& \mbox{minimize} & \frac{1}{2}\mathbf{x}^{T}P\mathbf{x} + \mathbf{q}^{T}\mathbf{x} + r \\
& \mbox{subject to} & \mathbf{l} \leq A\mathbf{x} \leq \mathbf{u}
\f}
</center>
in \f$ n \f$ real variables \f$ \mathbf{x} = (x_0, \ldots, x_{n-1}) \f$ and \f$ m \f$ constraints.
Here,
<UL>
<LI>\f$ P \f$ is a symmetric positive-semidefinite \f$ n \times n\f$ matrix (the quadratic objective function),
<LI>\f$ \mathbf{q} \f$ is an \f$ n \f$-dimensional vector (the linear objective function),
<LI>\f$ r \f$ is a constant,
<LI>\f$ A \f$ is an \f$ m\times n\f$ matrix (the constraint matrix),
<LI>\f$ \mathbf{l} \f$ is an \f$ m \f$-dimensional vector of lower constraint bounds,
where \f$ l_i \in \mathbb{R} \cup \{-\infty\} \f$ for all \f$ i \f$,
<LI>\f$ \mathbf{u} \f$ is an \f$ m \f$-dimensional vector of upper constraint bounds,
where \f$ u_i \in \mathbb{R} \cup \{+\infty\} \f$ for all \f$ i \f$.
</UL>
\cgalHasModel
`CGAL::OSQP_quadratic_program_traits<T>`
*/
class QuadraticProgramTraits {
public:
/// \name Memory
/// @{
/*!
Allocates memory for `n` variables and `m` constraints in `qp`.
*/
void resize(const std::size_t n, const std::size_t m) { }
/// @}
/// \name Initialization
/// @{
/*!
Sets the entries `Pij` and `Pji` of `qp` to `value`.
*/
void set_P(const std::size_t i, const std::size_t j, const FT value) { }
/*!
Sets the entry `qi` of `qp` to `value`.
*/
void set_q(const std::size_t i, const FT value) { }
/*!
Sets the entry `r` of `qp` to `value`.
*/
void set_r(const FT value) { }
/*!
Sets the entry `Aij` in the row `i` and column `j` of
the constraint matrix `A` of `qp` to `value`.
*/
void set_A(const std::size_t i, const std::size_t j, const FT value) { }
/*!
Sets the entry `li` of `qp` to `value`.
*/
void set_l(const std::size_t i, const FT value) { }
/*!
Sets the entry `ui` of `qp` to `value`.
*/
void set_u(const std::size_t i, const FT value) { }
/// @}
/// \name Solution
/// @{
/*!
\brief solves the quadratic program.
Number of values in `solution` equals to the number `n` of values in the vector `x`.
\tparam OutIterator
a model of `OutputIterator` that accepts values of type `FieldNumberType`
\param solution
an output iterator with the solution
\returns a status of the computation `success == true`
*/
template<typename OutIterator>
bool solve(OutIterator solution) { }
/// @}
};

53
Solver_interface/doc/Solver_interface/PackageDescription.txt
View File

@ -1,17 +1,35 @@
/// \defgroup PkgSolverInterfaceRef CGAL and Solvers Reference
/// \defgroup PkgSolverInterfaceConcepts Concepts
/// \ingroup PkgSolverInterfaceRef
///
/*!
\defgroup PkgSolverInterfaceRef CGAL and Solvers Reference
\defgroup PkgSolverInterfaceConcepts Concepts
\ingroup PkgSolverInterfaceRef
Concepts that describe various solvers.
\defgroup PkgSolverInterfaceLS Linear Systems
\ingroup PkgSolverInterfaceRef
Classes to define and solve linear systems.
\defgroup PkgSolverInterfaceMIP Mixed Integer Programming
\ingroup PkgSolverInterfaceRef
Classes to define and solve mixed integer programs.
\defgroup PkgSolverInterfaceNLP Nonlinear Programming
\ingroup PkgSolverInterfaceRef
Classes to define and solve nonlinear programs.
\addtogroup PkgSolverInterfaceRef
\cgalPkgDescriptionBegin{CGAL and Solvers,PkgSolverInterface}
\cgalPkgDescriptionBegin{CGAL and Solvers, PkgSolverInterface}
\cgalPkgPicture{solver.png}
\cgalPkgSummaryBegin
\cgalPkgAuthors{Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, Laurent Saboret, and Liangliang Nan}
\cgalPkgDesc{This package provides concepts and models for solving linear systems with dense or sparse matrices, Mixed Integer Programming (MIP) problems with or without constraints.}
\cgalPkgManuals{Chapter_CGAL_and_Solvers,PkgSolverInterfaceRef}
\cgalPkgDesc{This package provides concepts and models for solving linear systems with
dense or sparse matrices, mixed integer programming problems with or without constraints,
and nonlinear programming problems with or without constraints.}
\cgalPkgManuals{Chapter_CGAL_and_Solvers, PkgSolverInterfaceRef}
\cgalPkgSummaryEnd
\cgalPkgShortInfoBegin
\cgalPkgSince{4.8}
@ -23,16 +41,16 @@
\cgalClassifedRefPages
\cgalCRPSection{Concepts}
- `DiagonalizeTraits`
- `NormalEquationSparseLinearAlgebraTraits_d`
- `SparseLinearAlgebraTraits_d`
- `SparseLinearAlgebraWithFactorTraits_d`
- `SvdTraits`
- `MixedIntegerProgramTraits`
- `LinearProgramTraits`
- `QuadraticProgramTraits`
\cgalCRPSection{Classes}
\cgalCRPSection{Linear Systems}
- `CGAL::Eigen_solver_traits`
- `CGAL::Eigen_diagonalize_traits`
- `CGAL::Eigen_vector`
@ -40,10 +58,15 @@
- `CGAL::Eigen_sparse_matrix`
- `CGAL::Eigen_sparse_symmetric_matrix`
- `CGAL::Eigen_svd`
- `CGAL::Mixed_integer_program_traits`
- `CGAL::GLPK_mixed_integer_program_traits`
- `CGAL::SCIP_mixed_integer_program_traits`
\cgalCRPSection{Mixed Integer Programming}
- `CGAL::Variable`
- `CGAL::Linear_constraint`
- `CGAL::Linear_objective`
- `CGAL::Mixed_integer_program_traits`
- `CGAL::GLPK_mixed_integer_program_traits`
- `CGAL::SCIP_mixed_integer_program_traits`
\cgalCRPSection{Nonlinear Programming}
- `CGAL::OSQP_quadratic_program_traits`
*/

58
Solver_interface/doc/Solver_interface/Solver_interface.txt
View File

@ -7,9 +7,9 @@
\cgalAutoToc
\authors Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, Laurent Saboret, and Liangliang Nan
Several \cgal packages have to solve linear systems with dense or
sparse matrices, linear integer programs. This package provides
concepts and models for that purpose.
Several \cgal packages have to solve linear systems with dense or sparse matrices,
linear integer programs, and quadratic programs. This package provides concepts and
models for that purpose.
For linear systems, we generally provide models using the
\ref thirdpartyEigen library. Wrappers for the Eigen classes
@ -20,29 +20,27 @@ href="https://software.intel.com/en-us/mkl">Intel Math Kernel
Library (MKL)</a>.
For mixed integer programs (either constrained or unconstrained),
we provide models using \ref thirdpartySCIP and \ref thirdpartyGLPK
libraries. It is also possible to derive new models from other
high performance libraries, e.g.,
<a href = "https://projects.coin-or.org/Cbc"> CBC </a>,
<a href = "https://www.gurobi.com/"> Gurobi </a>.
we provide models using the \ref thirdpartySCIP and \ref thirdpartyGLPK
libraries.
For linear and quadratic programs, \cgal provides the built-in
\ref PkgQPSolver "CGAL Linear and Quadratic Programming Solver"
and we also provide a model using the \ref thirdpartyOSQP library.
\section SectionSolverDiagonalize Matrix Diagonalization
The concept `DiagonalizeTraits<T,dim>` defines an interface for the
The concept `DiagonalizeTraits<T, dim>` defines an interface for the
diagonalization and computation of eigenvectors and eigenvalues of a
symmetric matrix. `T` is the number type and `dim` is the dimension of
the matrices and vector (set to 3 by default). We provide the model
`Eigen_diagonalize_traits<T,dim>` that uses the \ref thirdpartyEigen library.
`Eigen_diagonalize_traits<T, dim>` that uses the \ref thirdpartyEigen library.
This is an example of an eigen decomposition of a matrix using this
class:
This is an example of an eigen decomposition of a matrix using this class:
\cgalExample{Solver_interface/diagonalize_matrix.cpp}
\section SectionSolverSVD Singular Value Decomposition
The concept `SvdTraits` defines an interface for solving in the least
@ -50,13 +48,12 @@ square sense a linear system with a singular value decomposition. The
field type is `double`. We provide the model `Eigen_svd` that uses the
\ref thirdpartyEigen library.
Here is a simple example that shows how to handle matrices, vectors
Here is a simple example that shows how to handle matrices, vectors,
and this solver:
\cgalExample{Solver_interface/singular_value_decomposition.cpp}
\section SectionSolverSparse Sparse Solvers
We define 3 concepts for sparse linear algebra:
@ -78,24 +75,22 @@ and solver is required:
typedef CGAL::Eigen_sparse_matrix<double>::EigenType EigenMatrix;
//iterative general solver
// iterative general solver
typedef CGAL::Eigen_solver_traits< Eigen::BiCGSTAB<EigenMatrix> > Iterative_general_solver;
//iterative symmetric solver
// iterative symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::ConjugateGradient<EigenMatrix> > Iterative_symmetric_solver;
//direct symmetric solver
// direct symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky<EigenMatrix> > Direct_symmetric_solver;
\endcode
Here is an example that shows how to fill the sparse matrix and call
the solver:
Here is an example that shows how to fill the sparse matrix and call the solver:
\cgalExample{Solver_interface/sparse_solvers.cpp}
\section SectionMIPSolver Mixed Integer Program Solvers
The concept `MixedIntegerProgramTraits` defines an interface for
@ -107,11 +102,23 @@ The field type is `double`. We provide two models of this concept:
`CGAL::SCIP_mixed_integer_program_traits` using \ref thirdpartySCIP.
Here is an example that shows how to formulate and solve a simple
mixed integer program using this solver.
mixed integer program using this solver:
\cgalExample{Solver_interface/mixed_integer_program.cpp}
\section SectionQPSolver Quadratic Program Solvers
The concept `QuadraticProgramTraits` defines an interface for quadratic programs (QP)
whereas a similar concept `LinearProgramTraits` gives an interface for linear programs (LP).
The model `CGAL::OSQP_quadratic_program_traits` provides a way to solve convex quadratic programs
with the dense or sparse interface.
Here is an example that shows how to formulate and solve a simple convex quadratic
program using the latter solver:
\cgalExample{Solver_interface/osqp_quadratic_program.cpp}
\section SolversHistory Implementation History
This package is the result of the increasing needs for linear solvers
@ -125,10 +132,11 @@ PkgSurfaceMeshSkeletonization and \ref PkgSurfaceMeshDeformation
extended the existing concepts. Liangliang Nan introduced the concept
`MixedIntegerProgramTraits` and two models for solving mixed integer
programs when implementing the \ref PkgPolygonalSurfaceReconstruction
package.
package. The concepts and models for solving linear and quadratic programs
were later introduced by Dmitry Anisimov and Mael Rouxel-Labbé.
Simon Giraudot was responsible for gathering all concepts and classes, and also wrote this user manual
with the help of Andreas Fabri.
Simon Giraudot was responsible for gathering all concepts and classes,
and also wrote this user manual with the help of Andreas Fabri.
*/
} /* namespace CGAL */

1
Solver_interface/doc/Solver_interface/dependencies
View File

@ -9,3 +9,4 @@ Surface_mesh_deformation
Jet_fitting_3
Poisson_surface_reconstruction_3
Polygonal_surface_reconstruction
QP_solver

1
Solver_interface/doc/Solver_interface/examples.txt
View File

@ -3,4 +3,5 @@
\example Solver_interface/singular_value_decomposition.cpp
\example Solver_interface/sparse_solvers.cpp
\example Solver_interface/mixed_integer_program.cpp
\example Solver_interface/osqp_quadratic_program.cpp
*/

17
Solver_interface/examples/Solver_interface/CMakeLists.txt
View File

@ -19,13 +19,27 @@ if(TARGET CGAL::Eigen3_support)
target_link_libraries(diagonalize_matrix PUBLIC CGAL::Eigen3_support)
endif()
create_single_source_cgal_program("mixed_integer_program.cpp")
find_package(OSQP QUIET)
include(CGAL_OSQP_support)
if(TARGET CGAL::OSQP_support)
create_single_source_cgal_program("osqp_quadratic_program.cpp")
target_link_libraries(osqp_quadratic_program PUBLIC CGAL::OSQP_support)
message("OSQP found and used")
else()
message(STATUS "NOTICE: OSQP was not found. OSQP examples won't be available.")
endif()
find_package(SCIP QUIET)
include(CGAL_SCIP_support)
if(TARGET CGAL::SCIP_support)
create_single_source_cgal_program("mixed_integer_program.cpp")
target_link_libraries(mixed_integer_program PUBLIC CGAL::SCIP_support)
message("SCIP found and used")
@ -36,6 +50,7 @@ else()
if(TARGET CGAL::GLPK_support)
create_single_source_cgal_program("mixed_integer_program.cpp")
target_link_libraries(mixed_integer_program PUBLIC CGAL::GLPK_support)
message("GLPK found and used")

15
Solver_interface/examples/Solver_interface/mixed_integer_program.cpp
View File

@ -32,9 +32,6 @@ typedef CGAL::GLPK_mixed_integer_program_traits<double> M
#endif
#if defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)
typedef typename MIP_Solver::Variable Variable;
typedef typename MIP_Solver::Linear_objective Linear_objective;
typedef typename MIP_Solver::Linear_constraint Linear_constraint;
@ -99,15 +96,3 @@ int main()
return EXIT_FAILURE;
}
}
#else
int main(int , char**)
{
std::cerr << "This test requires either GLPK or SCIP.\n";
return EXIT_SUCCESS;
}
#endif // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)

65
Solver_interface/examples/Solver_interface/osqp_quadratic_program.cpp Normal file
View File

@ -0,0 +1,65 @@
/*
* This example shows how to formulate and solve the following QP problem:
* https://osqp.org/docs/examples/setup-and-solve.html
*
* Minimize
* Objective:
* 1/2(4x1^2 + 2x1x2 + 2x2^2) + (x1 + x2) + 0 or in the matrix form
* 1/2 x^T|4 1|x + |1|^Tx + 0
* |1 2| |1|
* Subject to
* 1 <= x1 + x2 <= 1
* 0 <= x1 <= 0.7
* 0 <= x2 <= 0.7 or in the matrix form
* |1| |1 1| |1.0|
* |0| <= |1 0|x <= |0.7|
* |0| |0 1| |0.7|
*
* Expected results: x1 = 0.3; x2 = 0.7;
*/
#include <vector>
#include <iostream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/OSQP_quadratic_program_traits.h>
using Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
int main(void) {
const std::size_t n = 2; // number of variables
const std::size_t m = 3; // number of constraints
CGAL::OSQP_quadratic_program_traits<FT> osqp(n, m);
osqp.set_P(0, 0, 4);
osqp.set_P(0, 1, 1);
osqp.set_P(1, 1, 2);
osqp.set_q(0, 1);
osqp.set_q(1, 1);
osqp.set_r(0);
osqp.set_A(0, 0, 1);
osqp.set_A(0, 1, 1);
osqp.set_A(1, 0, 1);
osqp.set_A(2, 1, 1);
osqp.set_l(0, 1);
osqp.set_l(1, 0);
osqp.set_l(2, 0);
osqp.set_u(0, 1);
osqp.set_u(1, 0.7);
osqp.set_u(2, 0.7);
std::vector<FT> x; x.reserve(2);
osqp.solve(std::back_inserter(x));
std::cout << "solution (x1 x2): ";
for (const FT value : x) {
std::cout << value << "; ";
}
std::cout << std::endl;
}

2
Solver_interface/include/CGAL/Default_diagonalize_traits.h
View File

@ -21,7 +21,7 @@
namespace CGAL {
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceLS
///
/// The class `Default_diagonalize_traits` is a wrapper designed to automatically
/// use `Eigen_diagonalize_traits` if Eigen is available and otherwise use

2
Solver_interface/include/CGAL/Diagonalize_traits.h
View File

@ -29,7 +29,7 @@ lead to precision issues, please use CGAL::Eigen_diagonalize_traits")
namespace CGAL {
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceLS
///
/// The class `Diagonalize_traits` provides an internal
/// implementation for the diagonalization of Variance-Covariance

2
Solver_interface/include/CGAL/Eigen_diagonalize_traits.h
View File

@ -40,7 +40,7 @@ struct Restricted_FT<float> { typedef float type; };
}
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceLS
///
/// The class `Eigen_diagonalize_traits` provides an interface to the
/// diagonalization of covariance matrices of \ref thirdpartyEigen

2
Solver_interface/include/CGAL/Eigen_matrix.h
View File

@ -20,7 +20,7 @@
namespace CGAL {
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceLS
The class `Eigen_matrix` is a wrapper around `Eigen` matrix type
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html">`Eigen::Matrix`</a>.

6
Solver_interface/include/CGAL/Eigen_solver_traits.h
View File

@ -63,14 +63,14 @@ struct Get_eigen_matrix< ::Eigen::SparseLU<EigenMatrix, EigenOrdering >, FT>
} //internal
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceLS
The class `Eigen_solver_traits` provides an interface to the sparse solvers of \ref thirdpartyEigen "Eigen".
\ref thirdpartyEigen "Eigen" version 3.1 (or later) must be available on the system.
\cgalModels `SparseLinearAlgebraWithFactorTraits_d` and `NormalEquationSparseLinearAlgebraTraits_d`
\tparam EigenSolverT A sparse solver of \ref thirdpartyEigen "Eigen". The default solver is the iterative bi-congugate gradient stabilized solver `Eigen::BiCGSTAB` for `double`.
\tparam EigenSolverT A sparse solver of \ref thirdpartyEigen "Eigen". The default solver is the iterative bi-conjugate gradient stabilized solver `Eigen::BiCGSTAB` for `double`.
\sa `CGAL::Eigen_sparse_matrix<T>`
\sa `CGAL::Eigen_sparse_symmetric_matrix<T>`
@ -230,7 +230,7 @@ protected:
};
// Specialization of the solver for BiCGSTAB as for surface parameterization,
// the intializer should be a vector of one's (this was the case in 3.1-alpha
// the initializer should be a vector of one's (this was the case in 3.1-alpha
// but not in the official 3.1).
template<>
class Eigen_solver_traits<Eigen::BiCGSTAB<Eigen_sparse_matrix<double>::EigenType> >

4
Solver_interface/include/CGAL/Eigen_sparse_matrix.h
View File

@ -18,7 +18,7 @@
namespace CGAL {
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceLS
The class `Eigen_sparse_matrix` is a wrapper around `Eigen` matrix type
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix`</a>
@ -298,7 +298,7 @@ private:
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceRefLS
The class `Eigen_sparse_symmetric_matrix` is a wrapper around `Eigen` matrix type
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html">`Eigen::SparseMatrix` </a>

2
Solver_interface/include/CGAL/Eigen_svd.h
View File

@ -25,7 +25,7 @@
namespace CGAL {
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceLS
The class `Eigen_svd` provides an algorithm to solve in the least
square sense a linear system with a singular value decomposition using

2
Solver_interface/include/CGAL/Eigen_vector.h
View File

@ -17,7 +17,7 @@
namespace CGAL {
/*!
\ingroup PkgSolverInterfaceRef
\ingroup PkgSolverInterfaceLS
The class `Eigen_vector` is a wrapper around `Eigen`
<a href="http://eigen.tuxfamily.org/dox/classEigen_1_1Matrix.html">vector type</a>,

2
Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h
View File

@ -46,7 +46,7 @@ int bound_type(FT lb, FT ub)
} // namespace internal
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// This class provides an interface for formulating and solving
/// constrained or unconstrained mixed integer programs using

8
Solver_interface/include/CGAL/Mixed_integer_program_traits.h
View File

@ -96,7 +96,7 @@ namespace CGAL {
};
/// \endcond
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// The variable of a mixed integer program.
///
@ -195,7 +195,7 @@ namespace CGAL {
};
/// \endcond
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// The linear constraint of a mixed integer program.
///
@ -227,7 +227,7 @@ namespace CGAL {
};
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// The linear objective of a mixed integer program.
///
@ -263,7 +263,7 @@ namespace CGAL {
/// \endcond
};
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// The class `CGAL::Mixed_integer_program_traits` provides an interface for
/// formulating and solving (constrained or unconstrained) mixed integer

380
Solver_interface/include/CGAL/OSQP_quadratic_program_traits.h Normal file
View File

@ -0,0 +1,380 @@
// Copyright (c) 2020-2021 GeometryFactory SARL (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org)
//
// $URL$
// $Id$
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Dmitry Anisimov
// Mael Rouxel-Labbé
//
#ifndef CGAL_OSQP_QUADRATIC_PROGRAM_TRAITS_H
#define CGAL_OSQP_QUADRATIC_PROGRAM_TRAITS_H
#if defined(CGAL_USE_OSQP) || defined(DOXYGEN_RUNNING)
// STL includes.
#include <tuple>
#include <vector>
#include <utility>
#include <exception>
#include <memory>
// CGAL includes.
#include <CGAL/assertions.h>
// OSQP includes.
#include <osqp/osqp.h>
namespace CGAL {
/*!
\ingroup PkgSolverInterfaceNLP
\brief wraps the external OSQP solver.
This class provides an interface for formulating and solving
constrained or unconstrained quadratic programs using the \ref thirdpartyOSQP "OSQP"
library, which must be available on the system.
\tparam FT
number type that `FieldNumberType`
\note The `FT` type is provided for convenience. Internally, this FT type is converted
to `c_float` type that can be set either to `float` or `double`. By default, the `double`
type is used. After the optimization is complete, the `c_float` type is converted back to `FT`.
See more about `c_float` <a href="https://osqp.org/docs/interfaces/cc++#data-types">here</a>.
\cgalModels `QuadraticProgramTraits`
*/
template<typename FT>
class OSQP_quadratic_program_traits
{
using Triplet = std::tuple<std::size_t, std::size_t, FT>; // row, col, value
private:
std::size_t n; // number of variables
std::size_t m; // number of constraints
std::vector<Triplet> P_vec, A_vec;
std::vector<FT> q_vec, l_vec, u_vec;
public:
/// %Default constructor
OSQP_quadratic_program_traits() : n(0), m(0) { }
/// Constructor
/// \param n the number of variables
OSQP_quadratic_program_traits(const std::size_t n)
: n(n), m(0)
{
// 6*n, just guessing that this is some kind of sparse problem on a nice 2D mesh
P_vec.reserve(6 * n);
q_vec.reserve(n);
}
/// Constructor
/// \param n the number of variables
/// \param m the number of constraints
OSQP_quadratic_program_traits(const std::size_t n, const std::size_t m)
: n(n), m(m)
{
P_vec.reserve(6 * n);
q_vec.reserve(n);
A_vec.reserve(m);
l_vec.reserve(m);
u_vec.reserve(m);
}
public:
/// Resets the problem, removing all existing entries and setting sizes to `0`.
void clear()
{
n = m = 0;
P_vec.clear();
A_vec.clear();
q_vec.clear();
l_vec.clear();
u_vec.clear();
}
/// Changes the number of variables and the number of constraints of the problem.
///
/// \warning Calling this function also clears all previous entries.
void resize(const int new_n,
const int new_m = 0)
{
clear();
n = new_n;
m = new_m;
P_vec.reserve(6 * n);
q_vec.reserve(n);
if(m > 0)
{
A_vec.reserve(m);
l_vec.reserve(m);
u_vec.reserve(m);
}
}
/// \cond SKIP_IN_MANUAL
void set_P(const std::size_t i, const std::size_t j, const FT value)
{
if(j < i)
return;
if(j >= n) // no need to test i since j >= i
n = j+1;
P_vec.emplace_back(i, j, value);
}
void set_q(const std::size_t i, const FT value)
{
if(i >= n)
n = i+1;
if(i >= q_vec.size())
q_vec.resize(n);
q_vec[i] = value;
}
void set_r(const FT) {
// is not used here
}
void set_A(const std::size_t i, const std::size_t j, const FT value)
{
if(i >= m)
m = i+1;
if(j >= n)
n = j+1;
A_vec.emplace_back(i, j, value);
}
void set_l(const std::size_t i, const FT value)
{
if(i >= m)
m = i+1;
if(i >= l_vec.size())
l_vec.resize(m);
l_vec[i] = value;
}
void set_u(const std::size_t i, const FT value)
{
if(i >= m)
m = i+1;
if(i >= u_vec.size())
u_vec.resize(m);
u_vec[i] = value;
}
template<typename OutIterator>
bool solve(OutIterator solution,
const double eps_abs = 1e-10,
const double eps_rel = 1e-15,
const bool verbose = false)
{
if(verbose)
{
std::cout << "num variables = " << n << std::endl;
std::cout << "num constraints = " << m << std::endl;
}
CGAL_precondition(n >= 1); // m >= 0
CGAL_precondition(q_vec.size() == n);
CGAL_precondition(l_vec.size() == m && l_vec.size() == m);
const c_int P_nnz = static_cast<c_int>(P_vec.size());
auto P_x = std::make_unique<c_float[]>(P_nnz);
auto P_i = std::make_unique<c_int[]>(P_nnz);
auto P_p = std::make_unique<c_int[]>(n + 1);
set_matrix_from_triplets("P", P_vec, P_x.get(), P_i.get(), P_p.get());
if(verbose) std::cout << "P_nnz: " << P_nnz << std::endl;
const c_int A_nnz = static_cast<c_int>(A_vec.size());
auto A_x = std::make_unique<c_float[]>(A_nnz);
auto A_i = std::make_unique<c_int[]>(A_nnz);
auto A_p = std::make_unique<c_int[]>(n + 1);
set_matrix_from_triplets("A", A_vec, A_x.get(), A_i.get(), A_p.get());
if(verbose) std::cout << "A_nnz: " << A_nnz << std::endl;
const c_int q_size = static_cast<c_int>(q_vec.size());
const c_int l_size = static_cast<c_int>(l_vec.size());
const c_int u_size = static_cast<c_int>(u_vec.size());
auto q_x = std::make_unique<c_float[]>(q_size);
auto l_x = std::make_unique<c_float[]>(l_size);
auto u_x = std::make_unique<c_float[]>(u_size);
set_qlu_data(q_x.get(), l_x.get(), u_x.get());
// Problem settings.
OSQPSettings *settings = (OSQPSettings *) malloc(sizeof(OSQPSettings));
CGAL_assertion(settings);
// Structures.
OSQPWorkspace *work;
OSQPData *data = (OSQPData *) malloc(sizeof(OSQPData));
CGAL_assertion(data);
// Populate data.
data->n = static_cast<c_int>(n);
data->m = static_cast<c_int>(m);
data->P = csc_matrix(data->n, data->n, P_nnz, P_x.get(), P_i.get(), P_p.get());
CGAL_assertion(data->P);
data->q = q_x.get();
data->A = csc_matrix(data->m, data->n, A_nnz, A_x.get(), A_i.get(), A_p.get());
CGAL_assertion(data->A);
data->l = l_x.get();
data->u = u_x.get();
// Set solver settings.
osqp_set_default_settings(settings);
settings->eps_abs = eps_abs;
settings->eps_rel = eps_rel;
settings->verbose = verbose;
// Set workspace.
osqp_setup(&work, data, settings);
// Solve problem.
c_int exitflag = -1;
try
{
exitflag = osqp_solve(work);
}
catch(std::exception& e)
{
std::cerr << "ERROR: OSQP solver has thrown an exception!" << std::endl;
std::cerr << e.what() << std::endl;
}
const bool success = (exitflag == 0);
// Create solution.
const c_float *x = work->solution->x;
for(std::size_t i=0; i<n; ++i)
{
const FT value{x[i]};
*(++solution) = value;
}
// Clean workspace.
osqp_cleanup(work);
c_free(data->A);
c_free(data->P);
c_free(data);
c_free(settings);
return success;
}
/// \endcond
private:
// Based on the code in scipy, function coo_tocsr()
void set_matrix_from_triplets(const std::string /* name */,
const std::vector<Triplet>& triplets,
c_float *M_x,
c_int *M_i,
c_int *M_p) const
{
const std::size_t nnz = triplets.size();
// Compute the number of non-zero entries in each column of the sparse matrix A
std::fill(M_p, M_p + n, 0);
for(std::size_t k=0; k<nnz; ++k)
{
M_p[std::get<1>(triplets[k])]++;
}
// Fill M_p
c_int cumsum = 0;
for(std::size_t j=0; j<n; ++j)
{
const c_int tmp = M_p[j];
M_p[j] = cumsum;
cumsum += tmp;
}
M_p[n] = nnz;
// Write Ai, Ax into M_i, M_x
for(std::size_t k=0; k<nnz; ++k)
{
const c_int col = static_cast<c_int>(std::get<1>(triplets[k]));
const c_int dest = M_p[col];
M_i[dest] = static_cast<c_int>(std::get<0>(triplets[k]));
M_x[dest] = c_float(CGAL::to_double(std::get<2>(triplets[k])));
M_p[col]++;
}
c_int last = 0;
for(std::size_t j=0; j<=n; ++j)
{
const c_int tmp = M_p[j];
M_p[j] = last;
last = tmp;
}
// std::cout << name + "_x: ";
// for(std::size_t i=0; i<nnz; ++i)
// std::cout << M_x[i] << " ";
// std::cout << std::endl;
// std::cout << name + "_i: ";
// for(std::size_t i=0; i<nnz; ++i)
// std::cout << M_i[i] << " ";
// std::cout << std::endl;
// std::cout << name + "_p: ";
// for(std::size_t i=0; i<(n+1); ++i)
// std::cout << M_p[i] << " ";
// std::cout << std::endl;
}
void set_qlu_data(c_float *q_x,
c_float *l_x,
c_float *u_x) const
{
for(std::size_t i=0; i<n; ++i)
{
q_x[i] = c_float(CGAL::to_double(q_vec[i]));
}
for(std::size_t i=0; i<m; ++i)
{
l_x[i] = c_float(CGAL::to_double(l_vec[i]));
u_x[i] = c_float(CGAL::to_double(u_vec[i]));
}
// std::cout << "q_x: ";
// for(std::size_t i=0; i<n; ++i)
// std::cout << q_x[i] << " ";
// std::cout << std::endl;
// std::cout << "l_x: ";
// for(std::size_t i=0; i<m; ++i)
// std::cout << l_x[i] << " ";
// std::cout << std::endl;
// std::cout << "u_x: ";
// for(std::size_t i=0; i<m; ++i)
// std::cout << u_x[i] << " ";
// std::cout << std::endl;
}
};
} // namespace CGAL
#endif // CGAL_USE_OSQP or DOXYGEN_RUNNING
#endif // CGAL_OSQP_QUADRATIC_PROGRAM_TRAITS_H

2
Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h
View File

@ -25,7 +25,7 @@
namespace CGAL {
/// \ingroup PkgSolverInterfaceRef
/// \ingroup PkgSolverInterfaceMIP
///
/// This class provides an interface for formulating and solving
/// constrained or unconstrained mixed integer programs using

4
Solver_interface/package_info/Solver_interface/description.txt
View File

@ -1,2 +1,2 @@
Solver_interface: contains classes to interface third party linear algebra solvers (such as Eigen for example),
mixed integer program solvers (such as SCIP for example) to CGAL.
This package contains classes to interface third party linear algebra solvers (such as Eigen for example),
mixed integer program solvers (such as SCIP for example), and nonlinear program solvers (such as OSQP for example).