diff --git a/Documentation/doc/Documentation/Third_party.txt b/Documentation/doc/Documentation/Third_party.txt index 64d2da2e150..66454eb448b 100644 --- a/Documentation/doc/Documentation/Third_party.txt +++ b/Documentation/doc/Documentation/Third_party.txt @@ -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 `https://www.gnu.org/software/glpk/`. @@ -309,12 +306,18 @@ The \glpk web site is `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 `http://scip.zib.de/`. +\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 `https://osqp.org`. + */ diff --git a/Documentation/doc/resources/1.8.13/BaseDoxyfile.in b/Documentation/doc/resources/1.8.13/BaseDoxyfile.in index 90cb15aabb7..259b0709b4a 100644 --- a/Documentation/doc/resources/1.8.13/BaseDoxyfile.in +++ b/Documentation/doc/resources/1.8.13/BaseDoxyfile.in @@ -266,6 +266,7 @@ ALIASES = "cgal=%CGAL" \ "ceres=Ceres" \ "glpk=GLPK" \ "scip=SCIP" \ + "osqp=OSQP" \ "rs=RS" \ "rs3=RS3" \ "unix=Unix" \ diff --git a/Documentation/doc/resources/1.8.14/BaseDoxyfile.in b/Documentation/doc/resources/1.8.14/BaseDoxyfile.in index 96fd516d1ca..f0586530fc4 100644 --- a/Documentation/doc/resources/1.8.14/BaseDoxyfile.in +++ b/Documentation/doc/resources/1.8.14/BaseDoxyfile.in @@ -267,6 +267,7 @@ ALIASES = "cgal=%CGAL" \ "ceres=Ceres" \ "glpk=GLPK" \ "scip=SCIP" \ + "osqp=OSQP" \ "rs=RS" \ "rs3=RS3" \ "unix=Unix" \ diff --git a/Documentation/doc/resources/1.8.20/BaseDoxyfile.in b/Documentation/doc/resources/1.8.20/BaseDoxyfile.in index afb52c9c5af..7cdd82581dc 100644 --- a/Documentation/doc/resources/1.8.20/BaseDoxyfile.in +++ b/Documentation/doc/resources/1.8.20/BaseDoxyfile.in @@ -289,6 +289,7 @@ ALIASES = "cgal=%CGAL" \ "ceres=Ceres" \ "glpk=GLPK" \ "scip=SCIP" \ + "osqp=OSQP" \ "rs=RS" \ "rs3=RS3" \ "unix=Unix" \ diff --git a/Documentation/doc/resources/1.8.4/BaseDoxyfile.in b/Documentation/doc/resources/1.8.4/BaseDoxyfile.in index 9124909bf7c..fcb7026306b 100644 --- a/Documentation/doc/resources/1.8.4/BaseDoxyfile.in +++ b/Documentation/doc/resources/1.8.4/BaseDoxyfile.in @@ -235,6 +235,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" diff --git a/Installation/cmake/modules/CGAL_OSQP_support.cmake b/Installation/cmake/modules/CGAL_OSQP_support.cmake new file mode 100644 index 00000000000..6252d3aeeae --- /dev/null +++ b/Installation/cmake/modules/CGAL_OSQP_support.cmake @@ -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() diff --git a/Installation/cmake/modules/FindOSQP.cmake b/Installation/cmake/modules/FindOSQP.cmake new file mode 100644 index 00000000000..1fe8e048769 --- /dev/null +++ b/Installation/cmake/modules/FindOSQP.cmake @@ -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() diff --git a/Solver_interface/doc/Solver_interface/Concepts/LinearProgramTraits.h b/Solver_interface/doc/Solver_interface/Concepts/LinearProgramTraits.h new file mode 100644 index 00000000000..9d44b4a667f --- /dev/null +++ b/Solver_interface/doc/Solver_interface/Concepts/LinearProgramTraits.h @@ -0,0 +1,112 @@ +/*! +\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: +
+\f{eqnarray*}{ +& \mbox{minimize} & \mathbf{q}^{T}\mathbf{x} + r \\ +& \mbox{subject to} & \mathbf{l} \leq A\mathbf{x} \leq \mathbf{u} +\f} +
+in \f$ n \f$ real variables \f$ \mathbf{x} = (x_0, \ldots, x_{n-1}) \f$ and \f$ m \f$ constraints. + +Here, + +*/ +class LinearProgramTraits { + +public: + + /// \name Memory + /// @{ + + /*! + Allocates memory for `n` values in the vector `q`. + */ + void reserve_q(const std::size_t n) { } + + /*! + Allocates memory for `k` non-zero values in the matrix `A`. + */ + void reserve_A(const std::size_t k) { } + + /*! + Allocates memory for `m` values in the vector `l`. + */ + void reserve_l(const std::size_t m) { } + + /*! + Allocates memory for `m` values in the vector `u`. + */ + void reserve_u(const std::size_t m) { } + + /// @} + + /// \name Initialization + /// @{ + + /*! + Sets the entry `qj` of `lp` to `value`. + */ + void set_q(const std::size_t j, 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 `FT` + + \param solution + an output iterator with the solution + + \returns a status of the computation `success == true` + */ + template + bool solve(OutIterator solution) + { } + + /// @} +}; diff --git a/Solver_interface/doc/Solver_interface/Concepts/QuadraticProgramTraits.h b/Solver_interface/doc/Solver_interface/Concepts/QuadraticProgramTraits.h new file mode 100644 index 00000000000..504bdcb1296 --- /dev/null +++ b/Solver_interface/doc/Solver_interface/Concepts/QuadraticProgramTraits.h @@ -0,0 +1,129 @@ +/*! +\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: +
+\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} +
+in \f$ n \f$ real variables \f$ \mathbf{x} = (x_0, \ldots, x_{n-1}) \f$ and \f$ m \f$ constraints. + +Here, + + +\cgalHasModel +`CGAL::OSQP_quadratic_program_traits` +*/ +class QuadraticProgramTraits { + +public: + + /// \name Memory + /// @{ + + /*! + Allocates memory for `k` non-zero values in the matrix `P`. + */ + void reserve_P(const std::size_t k) { } + + /*! + Allocates memory for `n` values in the vector `q`. + */ + void reserve_q(const std::size_t n) { } + + /*! + Allocates memory for `k` non-zero values in the matrix `A`. + */ + void reserve_A(const std::size_t k) { } + + /*! + Allocates memory for `m` values in the vector `l`. + */ + void reserve_l(const std::size_t m) { } + + /*! + Allocates memory for `m` values in the vector `u`. + */ + void reserve_u(const std::size_t m) { } + + /// @} + + /// \name Initialization + /// @{ + + /*! + Sets the entries `Pij` and `Pji` of `qp` to `value`. + + Note that you should define only the upper triangular part of the matrix! + */ + void set_P(const std::size_t i, const std::size_t j, const FT value) + { } + + /*! + Sets the entry `qj` of `qp` to `value`. + */ + void set_q(const std::size_t j, 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 `FT` + + \param solution + an output iterator with the solution + + \returns a status of the computation `success == true` + */ + template + bool solve(OutIterator solution) + { } + + /// @} +}; diff --git a/Solver_interface/doc/Solver_interface/PackageDescription.txt b/Solver_interface/doc/Solver_interface/PackageDescription.txt index dcc72642a69..55f61bed124 100644 --- a/Solver_interface/doc/Solver_interface/PackageDescription.txt +++ b/Solver_interface/doc/Solver_interface/PackageDescription.txt @@ -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` */ diff --git a/Solver_interface/doc/Solver_interface/Solver_interface.txt b/Solver_interface/doc/Solver_interface/Solver_interface.txt index e5fded45ffa..4a42e373609 100644 --- a/Solver_interface/doc/Solver_interface/Solver_interface.txt +++ b/Solver_interface/doc/Solver_interface/Solver_interface.txt @@ -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,31 @@ href="https://software.intel.com/en-us/mkl">Intel Math Kernel Library (MKL). For mixed integer programs (either constrained or unconstrained), -we provide models using \ref thirdpartySCIP and \ref thirdpartyGLPK +we provide models using the \ref thirdpartySCIP and \ref thirdpartyGLPK libraries. It is also possible to derive new models from other -high performance libraries, e.g., - CBC , - Gurobi . +high performance libraries such as CBC or +Gurobi. +For linear and quadratic programs, we provide a model using the built-in +\ref PkgQPSolver "CGAL Linear and Quadratic Programming Solver" and a model using +the external \ref thirdpartyOSQP library. It is also possible to derive new models +from other quadratic programming libraries such as +OOQP for example. \section SectionSolverDiagonalize Matrix Diagonalization -The concept `DiagonalizeTraits` defines an interface for the +The concept `DiagonalizeTraits` 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` that uses the \ref thirdpartyEigen library. +`Eigen_diagonalize_traits` 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 +52,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 +79,22 @@ and solver is required: typedef CGAL::Eigen_sparse_matrix::EigenType EigenMatrix; -//iterative general solver +// iterative general solver typedef CGAL::Eigen_solver_traits< Eigen::BiCGSTAB > Iterative_general_solver; -//iterative symmetric solver +// iterative symmetric solver typedef CGAL::Eigen_solver_traits< Eigen::ConjugateGradient > Iterative_symmetric_solver; -//direct symmetric solver +// direct symmetric solver typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky > 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 +106,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 +136,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. -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 */ diff --git a/Solver_interface/doc/Solver_interface/dependencies b/Solver_interface/doc/Solver_interface/dependencies index 2dd25947f0a..b00ac4fc10c 100644 --- a/Solver_interface/doc/Solver_interface/dependencies +++ b/Solver_interface/doc/Solver_interface/dependencies @@ -9,3 +9,4 @@ Surface_mesh_deformation Jet_fitting_3 Poisson_surface_reconstruction_3 Polygonal_surface_reconstruction +QP_solver diff --git a/Solver_interface/doc/Solver_interface/examples.txt b/Solver_interface/doc/Solver_interface/examples.txt index 680c5416f6b..64cf531f9e7 100644 --- a/Solver_interface/doc/Solver_interface/examples.txt +++ b/Solver_interface/doc/Solver_interface/examples.txt @@ -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 */ diff --git a/Solver_interface/examples/Solver_interface/CMakeLists.txt b/Solver_interface/examples/Solver_interface/CMakeLists.txt index 8b1a72ef9b0..b700decef6d 100644 --- a/Solver_interface/examples/Solver_interface/CMakeLists.txt +++ b/Solver_interface/examples/Solver_interface/CMakeLists.txt @@ -20,6 +20,21 @@ if(TARGET CGAL::Eigen3_support) endif() create_single_source_cgal_program("mixed_integer_program.cpp") +create_single_source_cgal_program("osqp_quadratic_program.cpp") + +find_package(OSQP QUIET) +include(CGAL_OSQP_support) + +if(TARGET CGAL::OSQP_support) + + target_link_libraries(osqp_quadratic_program PUBLIC CGAL::OSQP_support) + message("OSQP found and used") + +else() + + message(STATUS "NOTICE: OSQP was not found. Complete OSQP examples won't be available.") + +endif() find_package(SCIP QUIET) include(CGAL_SCIP_support) diff --git a/Solver_interface/examples/Solver_interface/osqp_quadratic_program.cpp b/Solver_interface/examples/Solver_interface/osqp_quadratic_program.cpp new file mode 100644 index 00000000000..a7f72d8b3ee --- /dev/null +++ b/Solver_interface/examples/Solver_interface/osqp_quadratic_program.cpp @@ -0,0 +1,80 @@ +/* +* 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 +#include +#include + +#if defined(CGAL_USE_OSQP) +#include +#endif + +using Kernel = CGAL::Simple_cartesian; +using FT = typename Kernel::FT; + +#if defined(CGAL_USE_OSQP) + +int main(void) { + + CGAL::OSQP_quadratic_program_traits osqp; + + osqp.reserve_P(3); + osqp.set_P(0, 0, 4); + osqp.set_P(0, 1, 1); + osqp.set_P(1, 1, 2); + + osqp.reserve_q(2); + osqp.set_q(0, 1); + osqp.set_q(1, 1); + + osqp.set_r(0); + + osqp.reserve_A(4); + 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.reserve_l(3); + osqp.set_l(0, 1); + osqp.set_l(1, 0); + osqp.set_l(2, 0); + + osqp.reserve_u(3); + osqp.set_u(0, 1); + osqp.set_u(1, 0.7); + osqp.set_u(2, 0.7); + + std::vector 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; +} + +#else +int main(void) { + std::cout << "This example requires the OSQP library." << std::endl; + return EXIT_SUCCESS; +} +#endif // defined(CGAL_USE_OSQP) diff --git a/Solver_interface/include/CGAL/Default_diagonalize_traits.h b/Solver_interface/include/CGAL/Default_diagonalize_traits.h index cbf009379a8..930dd681ce2 100644 --- a/Solver_interface/include/CGAL/Default_diagonalize_traits.h +++ b/Solver_interface/include/CGAL/Default_diagonalize_traits.h @@ -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 diff --git a/Solver_interface/include/CGAL/Diagonalize_traits.h b/Solver_interface/include/CGAL/Diagonalize_traits.h index 929752edba0..c4754bed5e9 100644 --- a/Solver_interface/include/CGAL/Diagonalize_traits.h +++ b/Solver_interface/include/CGAL/Diagonalize_traits.h @@ -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 diff --git a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h index 38ee3c7c4c7..e0321f57395 100644 --- a/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h +++ b/Solver_interface/include/CGAL/Eigen_diagonalize_traits.h @@ -40,7 +40,7 @@ struct Restricted_FT { typedef float type; }; } -/// \ingroup PkgSolverInterfaceRef +/// \ingroup PkgSolverInterfaceLS /// /// The class `Eigen_diagonalize_traits` provides an interface to the /// diagonalization of covariance matrices of \ref thirdpartyEigen diff --git a/Solver_interface/include/CGAL/Eigen_matrix.h b/Solver_interface/include/CGAL/Eigen_matrix.h index 644ae53babb..4a34a022d0d 100644 --- a/Solver_interface/include/CGAL/Eigen_matrix.h +++ b/Solver_interface/include/CGAL/Eigen_matrix.h @@ -20,7 +20,7 @@ namespace CGAL { /*! -\ingroup PkgSolverInterfaceRef +\ingroup PkgSolverInterfaceLS The class `Eigen_matrix` is a wrapper around `Eigen` matrix type `Eigen::Matrix`. diff --git a/Solver_interface/include/CGAL/Eigen_solver_traits.h b/Solver_interface/include/CGAL/Eigen_solver_traits.h index c7efdfda81d..e2fc4f64c0d 100644 --- a/Solver_interface/include/CGAL/Eigen_solver_traits.h +++ b/Solver_interface/include/CGAL/Eigen_solver_traits.h @@ -63,7 +63,7 @@ struct Get_eigen_matrix< ::Eigen::SparseLU, 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. diff --git a/Solver_interface/include/CGAL/Eigen_sparse_matrix.h b/Solver_interface/include/CGAL/Eigen_sparse_matrix.h index 54d8663d5ca..ff83d2a65fd 100644 --- a/Solver_interface/include/CGAL/Eigen_sparse_matrix.h +++ b/Solver_interface/include/CGAL/Eigen_sparse_matrix.h @@ -18,7 +18,7 @@ namespace CGAL { /*! -\ingroup PkgSolverInterfaceRef +\ingroup PkgSolverInterfaceLS The class `Eigen_sparse_matrix` is a wrapper around `Eigen` matrix type `Eigen::SparseMatrix` @@ -298,7 +298,7 @@ private: /*! -\ingroup PkgSolverInterfaceRef +\ingroup PkgSolverInterfaceRefLS The class `Eigen_sparse_symmetric_matrix` is a wrapper around `Eigen` matrix type `Eigen::SparseMatrix` diff --git a/Solver_interface/include/CGAL/Eigen_svd.h b/Solver_interface/include/CGAL/Eigen_svd.h index 9697401afe8..0841976072a 100644 --- a/Solver_interface/include/CGAL/Eigen_svd.h +++ b/Solver_interface/include/CGAL/Eigen_svd.h @@ -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 diff --git a/Solver_interface/include/CGAL/Eigen_vector.h b/Solver_interface/include/CGAL/Eigen_vector.h index d113f85e78c..d080777a2f8 100644 --- a/Solver_interface/include/CGAL/Eigen_vector.h +++ b/Solver_interface/include/CGAL/Eigen_vector.h @@ -17,7 +17,7 @@ namespace CGAL { /*! -\ingroup PkgSolverInterfaceRef +\ingroup PkgSolverInterfaceLS The class `Eigen_vector` is a wrapper around `Eigen` vector type, diff --git a/Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h b/Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h index 8b6a8e28f65..1218679b04a 100644 --- a/Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h +++ b/Solver_interface/include/CGAL/GLPK_mixed_integer_program_traits.h @@ -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 diff --git a/Solver_interface/include/CGAL/Mixed_integer_program_traits.h b/Solver_interface/include/CGAL/Mixed_integer_program_traits.h index df11eea5377..1b495e95824 100644 --- a/Solver_interface/include/CGAL/Mixed_integer_program_traits.h +++ b/Solver_interface/include/CGAL/Mixed_integer_program_traits.h @@ -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 diff --git a/Solver_interface/include/CGAL/OSQP_quadratic_program_traits.h b/Solver_interface/include/CGAL/OSQP_quadratic_program_traits.h new file mode 100644 index 00000000000..31dbcb85744 --- /dev/null +++ b/Solver_interface/include/CGAL/OSQP_quadratic_program_traits.h @@ -0,0 +1,281 @@ +// Copyright (c) 2020 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 +// + +#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 +#include +#include + +// CGAL includes. +#include + +// OSQP includes. +#include + +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 + + \cgalModels `QuadraticProgramTraits` + */ + template + class OSQP_quadratic_program_traits { + // row, col, value + using Triplet = std::tuple; + + public: + /// \cond SKIP_IN_MANUAL + void reserve_P(const std::size_t k) { + P_vec.reserve(k); + } + + void reserve_q(const std::size_t n) { + q_vec.reserve(n); + } + + void reserve_A(const std::size_t k) { + A_vec.reserve(k); + } + + void reserve_l(const std::size_t m) { + l_vec.reserve(m); + } + + void reserve_u(const std::size_t m) { + u_vec.reserve(m); + } + + void set_P(const std::size_t i, const std::size_t j, const FT value) { + P_vec.push_back(std::make_tuple(i, j, value)); + } + + void set_q(const std::size_t, const FT value) { + q_vec.push_back(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) { + A_vec.push_back(std::make_tuple(i, j, value)); + } + + void set_l(const std::size_t, const FT value) { + l_vec.push_back(value); + } + + void set_u(const std::size_t, const FT value) { + u_vec.push_back(value); + } + + std::size_t P_size() const { + return P_vec.size(); + } + + std::size_t q_size() const { + return q_vec.size(); + } + + std::size_t A_size() const { + return A_vec.size(); + } + + std::size_t l_size() const { + return l_vec.size(); + } + + std::size_t u_size() const { + return u_vec.size(); + } + + template + bool solve(OutIterator solution) { + + const std::size_t num_cols = std::get<1>(P_vec.back()) + 1; + CGAL_precondition(q_vec.size() == num_cols); + CGAL_precondition(l_vec.size() == u_vec.size()); + + const c_int P_nnz = static_cast(P_vec.size()); + c_float P_x[P_nnz]; + c_int P_i[P_nnz]; + c_int P_p[P_nnz + 1]; + set_matrix_from_triplets("P", P_vec, P_x, P_i, P_p); + // std::cout << "P_nnz: " << P_nnz << std::endl; + + const c_int A_nnz = static_cast(A_vec.size()); + c_float A_x[A_nnz]; + c_int A_i[A_nnz]; + c_int A_p[P_nnz + 1]; + set_matrix_from_triplets("A", A_vec, A_x, A_i, A_p); + // std::cout << "A_nnz: " << A_nnz << std::endl; + + const c_int q_size = static_cast(q_vec.size()); + const c_int l_size = static_cast(l_vec.size()); + const c_int u_size = static_cast(u_vec.size()); + + c_float q_x[q_size]; + c_float l_x[l_size]; + c_float u_x[u_size]; + set_qlu_data(q_x, l_x, u_x); + + // Problem settings. + OSQPSettings *settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings)); + + // Structures. + OSQPWorkspace *work; + OSQPData *data; + + // Populate data. + const c_int n = q_size; // number of variables + const c_int m = l_size; // number of constraints + // std::cout << "n: " << n << std::endl; + // std::cout << "m: " << m << std::endl; + + data = (OSQPData *)c_malloc(sizeof(OSQPData)); + data->n = n; + data->m = m; + data->P = csc_matrix(n, n, P_nnz, P_x, P_i, P_p); + data->q = q_x; + data->A = csc_matrix(m, n, A_nnz, A_x, A_i, A_p); + data->l = l_x; + data->u = u_x; + + // Set solver settings. + osqp_set_default_settings(settings); + settings->eps_rel = 1.0e-15; + settings->verbose = false; + + // Set workspace. + osqp_setup(&work, data, settings); + + // Solve problem. + const int exitflag = osqp_solve(work); + const bool success = exitflag == 0 ? true : false; + + // Create solution. + const c_float *x = work->solution->x; + for (c_int i = 0; i < n; ++i) { + const FT value = static_cast(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: + std::vector P_vec, A_vec; + std::vector q_vec, l_vec, u_vec; + + void set_matrix_from_triplets( + const std::string /* name */, + const std::vector& triplets, + c_float *M_x, c_int *M_i, c_int *M_p) const { + + M_p[0] = 0; + std::size_t count = 0; + const std::size_t num_cols = std::get<1>(triplets.back()); + std::size_t ref_col = 0; + for (ref_col = 0; ref_col <= num_cols; ++ref_col) { + std::size_t num_rows = 0; + for (std::size_t i = 0; i < triplets.size(); ++i) { + const std::size_t row = std::get<0>(triplets[i]); + const std::size_t col = std::get<1>(triplets[i]); + const double value = CGAL::to_double(std::get<2>(triplets[i])); + + if (col == ref_col) { + M_i[count] = row; M_x[count] = value; + ++count; ++num_rows; + } + } + M_p[ref_col + 1] = M_p[ref_col] + num_rows; + } + + // std::cout << name + "_x: "; + // for (std::size_t i = 0; i < count; ++i) + // std::cout << M_x[i] << " "; + // std::cout << std::endl; + + // std::cout << name + "_i: "; + // for (std::size_t i = 0; i < count; ++i) + // std::cout << M_i[i] << " "; + // std::cout << std::endl; + + // std::cout << name + "_p: "; + // for (std::size_t i = 0; i <= ref_col; ++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 { + + CGAL_assertion(q_vec.size() > 1); + CGAL_assertion(l_vec.size() == u_vec.size()); + const std::size_t n = q_vec.size(); // number of variables + const std::size_t m = l_vec.size(); // number of constraints + + CGAL_assertion(n > 1); + for (std::size_t i = 0; i < n; ++i) + q_x[i] = CGAL::to_double(q_vec[i]); + + CGAL_assertion(m >= 0); + for (std::size_t i = 0; i < m; ++i) { + l_x[i] = CGAL::to_double(l_vec[i]); + u_x[i] = 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 diff --git a/Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h b/Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h index ecc9320e103..2b1d7927447 100644 --- a/Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h +++ b/Solver_interface/include/CGAL/SCIP_mixed_integer_program_traits.h @@ -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 diff --git a/Solver_interface/package_info/Solver_interface/description.txt b/Solver_interface/package_info/Solver_interface/description.txt index 90a6f6a11f7..f31a09e2cf3 100644 --- a/Solver_interface/package_info/Solver_interface/description.txt +++ b/Solver_interface/package_info/Solver_interface/description.txt @@ -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).