1 // Copyright (c) 1997-2007 ETH Zurich (Switzerland).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/QP_solver/include/CGAL/QP_solver/QP_functions_impl.h $
7 // $Id: QP_functions_impl.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Bernd Gaertner <gaertner@inf.ethz.ch>
12 #ifndef CGAL_QP_FUNCTIONS_IMPL_H
13 #define CGAL_QP_FUNCTIONS_IMPL_H
14
15 #include <CGAL/license/QP_solver.h>
16
17
18 #include <iostream>
19 #include <fstream>
20 #include <CGAL/iterator.h>
21 #include <CGAL/QP_solver/QP_solver.h>
22 #include <CGAL/QP_models.h>
23 #include <CGAL/QP_solution.h>
24
25 namespace CGAL {
26
27 namespace QP_functions_detail {
28 // test whether the system is of the form A x == b (equations only)
29 template <typename R>
is_in_equational_form(const R & r)30 bool is_in_equational_form (const R& r)
31 {
32 typename R::R_iterator it = r.get_r();
33 typename R::R_iterator end = it + r.get_m();
34 for (; it < end; ++it)
35 if (*it != CGAL::EQUAL) return false;
36 return true;
37 }
38
39 // test whether the row vectors of A that correpsond to equations
40 // are linearly independent; this is done using type ET. The value
41 // type of LinearInequalitySystem must be convertible to ET
42 template <class Ar, class ET>
has_linearly_independent_equations(const Ar & ar,const ET &)43 bool has_linearly_independent_equations
44 (const Ar& ar, const ET& /*dummy*/) {
45 // we solve the following auxiliary LP, using exact type ET:
46 // --------
47 // min 0
48 // A x r 0
49 // x >= 0
50 // --------
51 // Then A has linearly independent equations if and only if all
52 // artificials have left the basis after phase I; the QP_solver
53 // diagnostics tells us this
54 //
55 // auxiliary LP type
56 typedef typename
57 std::iterator_traits<typename Ar::C_iterator>::value_type C_value;
58 typedef typename
59 std::iterator_traits<typename Ar::B_iterator>::value_type B_value;
60 typedef Const_oneset_iterator <C_value> C_iterator;
61 typedef Const_oneset_iterator <B_value> B_iterator;
62 typedef Nonnegative_linear_program_from_iterators
63 <typename Ar::A_iterator, B_iterator,
64 typename Ar::R_iterator, C_iterator> LP;
65
66 // auxiliary LP
67 LP lp (ar.get_n(), ar.get_m(), ar.get_a(), B_iterator(0), ar.get_r(), C_iterator(0));
68
69 // solver Tags
70 typedef QP_solver_impl::QP_tags<
71 Tag_true, // Is_linear
72 Tag_true> // Is_nonnegative
73 Tags;
74
75 // solver type
76 typedef QP_solver<LP, ET, Tags> Solver;
77
78 // now solve auxiliary LP and compute predicate value
79 Solver solver (lp);
80 return !solver.diagnostics.redundant_equations;
81 }
82
83 // helper for MPS output: BOUNDS
84 template <typename P>
print_bounds(std::ostream &,const P &,CGAL::Tag_true)85 void print_bounds
86 (std::ostream& , const P& ,
87 CGAL::Tag_true /*is_nonnegative*/)
88 {
89 // nop (default bounds are nonnegative)
90 }
91
92 // helper for MPS output: BOUNDS
93 template <typename P>
print_bounds(std::ostream & out,const P & p,CGAL::Tag_false)94 void print_bounds
95 (std::ostream& out, const P& p,
96 CGAL::Tag_false /*is_nonnegative*/)
97 {
98 typename P::FL_iterator fl = p.get_fl();
99 typename P::FU_iterator fu = p.get_fu();
100 typename P::L_iterator l = p.get_l();
101 typename P::U_iterator u = p.get_u();
102 int n = p.get_n();
103 out << "BOUNDS\n";
104 for (int j=0; j<n; ++j, ++fl, ++l, ++fu, ++u) {
105 if (!*fl || !CGAL::is_zero(*l)) {
106 if (*fl)
107 out << " LO BND x" << j << " " << *l << "\n";
108 else
109 out << " MI BND x" << j << "\n";
110 }
111 if (*fu)
112 out << " UP BND x" << j << " " << *u << "\n";
113 }
114 }
115
116 // helper for MPS output: DMATRIX/QMATRIX
117 template <typename P>
print_qmatrix(std::ostream &,const P &,CGAL::Tag_true)118 void print_qmatrix
119 (std::ostream& , const P& ,
120 CGAL::Tag_true /*is_linear*/)
121 {
122 // nop
123 }
124
125 // helper for MPS output: DMATRIX/QMATRIX
126 template <typename P>
print_qmatrix(std::ostream & out,const P & p,CGAL::Tag_false)127 void print_qmatrix
128 (std::ostream& out, const P& p,
129 CGAL::Tag_false /*is_linear*/)
130 {
131 typename P::D_iterator it = p.get_d();
132 int n = p.get_n();
133 bool empty_D = true;
134 for (int i=0; i<n; ++i, ++it) {
135 // handle only entries on/below diagonal
136 for (int j=0; j<i+1; ++j)
137 if (!CGAL::is_zero(*(*it + j))) {
138 if (empty_D) {
139 // first time we see a nonzero entry
140 out << "QMATRIX\n";
141 empty_D = false;
142 }
143 out << " x" << i << " x" << j << " " << *(*it + j) << "\n";
144 // QMATRIX format prescribes symmetric matrix
145 if (i != j)
146 out << " x" << j << " x" << i << " " << *(*it + j) << "\n";
147 }
148 }
149 }
150
151 // check whether the two qp's have the same data; this is the case iff
152 // they agree in n, m, a, b, r, fl, l, fu, u, d, c, c0
153 // PRE: qp1, qp2 have the same internal number type
154 template <typename Quadratic_program1, typename Quadratic_program2>
are_equal_qp(const Quadratic_program1 & qp1,const Quadratic_program2 & qp2)155 bool are_equal_qp
156 (const Quadratic_program1 &qp1, const Quadratic_program2 &qp2)
157 {
158 bool return_val = true;
159 // check n
160 if (qp1.get_n() != qp2.get_n()) {
161 std::cerr << "Equality test fails with n: "
162 << qp1.get_n() << " vs. " << qp2.get_n() << std::endl;
163 return false; // wildly wrong, abort now
164 }
165 // check m
166 if (qp1.get_m() != qp2.get_m()) {
167 std::cerr << "Equality test fails with m: "
168 << qp1.get_m() << " vs. " << qp2.get_m() << std::endl;
169 return false; // wildly wrong, abort now
170 }
171 int n = qp1.get_n();
172 int m = qp1.get_m();
173 // check A
174 typename Quadratic_program1::A_iterator a1 = qp1.get_a();
175 typename Quadratic_program2::A_iterator a2 = qp2.get_a();
176 for (int j=0; j<n; ++j, ++a1, ++a2)
177 for (int i=0; i<m; ++i)
178 if (*((*a1)+i) != *((*a2)+i)) {
179 std::cerr << "Equality test fails with A["
180 << j << "][" << i << "]: "
181 << *((*a1)+i) << " vs. " << *((*a2)+i) << std::endl;
182 return_val = false;
183 }
184 // check b
185 typename Quadratic_program1::B_iterator b1 = qp1.get_b();
186 typename Quadratic_program2::B_iterator b2 = qp2.get_b();
187 for (int i=0; i<m; ++i, ++b1, ++b2)
188 if (*b1 != *b2) {
189 std::cerr << "Equality test fails with b[" << i << "]: "
190 << *b1 << " vs. " << *b2 << std::endl;
191 return_val = false;
192 }
193 // check r
194 typename Quadratic_program1::R_iterator r1 = qp1.get_r();
195 typename Quadratic_program2::R_iterator r2 = qp2.get_r();
196 for (int i=0; i<m; ++i, ++r1, ++r2)
197 if (*r1 != *r2) {
198 std::cerr << "Equality test fails with r[" << i << "]: "
199 << *r1 << " vs. " << *r2 << std::endl;
200 return_val = false;
201 }
202 // check fl, l
203 typename Quadratic_program1::FL_iterator fl1 = qp1.get_fl();
204 typename Quadratic_program2::FL_iterator fl2 = qp2.get_fl();
205 typename Quadratic_program1::L_iterator l1 = qp1.get_l();
206 typename Quadratic_program2::L_iterator l2 = qp2.get_l();
207 for (int j=0; j<n; ++j, ++fl1, ++fl2, ++l1, ++l2) {
208 if (*fl1 != *fl2) {
209 std::cerr << "Equality test fails with fl[" << j << "]: "
210 << *fl1 << " vs. " << *fl2 << std::endl;
211 return_val = false;
212 }
213 if ((*fl1 == true) && (*l1 != *l2)) {
214 std::cerr << "Equality test fails with l[" << j << "]: "
215 << *l1 << " vs. " << *l2 << std::endl;
216 return_val = false;
217 }
218 }
219
220 // check fu, u
221 typename Quadratic_program1::FU_iterator fu1 = qp1.get_fu();
222 typename Quadratic_program2::FU_iterator fu2 = qp2.get_fu();
223 typename Quadratic_program1::U_iterator u1 = qp1.get_u();
224 typename Quadratic_program2::U_iterator u2 = qp2.get_u();
225 for (int j=0; j<n; ++j, ++fu1, ++fu2, ++u1, ++u2) {
226 if (*fu1 != *fu2) {
227 std::cerr << "Equality test fails with fu[" << j << "]: "
228 << *fu1 << " vs. " << *fu2 << std::endl;
229 return_val = false;
230 }
231 if ((*fu1 == true) && (*u1 != *u2)) {
232 std::cerr << "Equality test fails with u[" << j << "]: "
233 << *u1 << " vs. " << *u2 << std::endl;
234 return_val = false;
235 }
236 }
237 // check d
238 typename Quadratic_program1::D_iterator d1 = qp1.get_d();
239 typename Quadratic_program2::D_iterator d2 = qp2.get_d();
240 for (int i=0; i<n; ++i, ++d1, ++d2)
241 for (int j=0; j<i+1; ++j) // only access entries on/below diagonal
242 if (*((*d1)+j) != *((*d2)+j)) {
243 std::cerr << "Equality test fails with D["
244 << i << "][" << j << "]: "
245 << *((*d1)+j) << " vs. " << *((*d2)+j) << std::endl;
246 return_val = false;
247 }
248 // check c
249 typename Quadratic_program1::C_iterator c1 = qp1.get_c();
250 typename Quadratic_program2::C_iterator c2 = qp2.get_c();
251 for (int j=0; j<n; ++j, ++c1, ++c2)
252 if (*c1 != *c2) {
253 std::cerr << "Equality test fails with c[" << j << "]: "
254 << *c1 << " vs. " << *c2 << std::endl;
255 return_val = false;
256 }
257 // check c0
258 typename Quadratic_program1::C_entry c01 = qp1.get_c0();
259 typename Quadratic_program2::C_entry c02 = qp2.get_c0();
260 if (c01 != c02) {
261 std::cerr << "Equality test fails with c0: "
262 << c01 << " vs. " << c02 << std::endl;
263 return_val = false;
264 }
265 return return_val;
266 }
267
268 template <typename P, typename Is_linear, typename Is_nonnegative>
print_program(std::ostream & out,const P & p,const std::string & problem_name,Is_linear is_linear,Is_nonnegative is_nonnegative)269 void print_program
270 (std::ostream& out, const P& p,
271 const std::string& problem_name,
272 Is_linear is_linear,
273 Is_nonnegative is_nonnegative)
274 {
275 // NAME:
276 out << "NAME " << problem_name << "\n";
277
278 int n = p.get_n();
279 int m = p.get_m();
280
281 // ROWS section:
282 typename P::R_iterator r = p.get_r();
283 out << "ROWS\n"
284 << " N obj\n"; // for the objective function
285 for (int i=0; i<m; ++i, ++r) {
286 if (*r == CGAL::SMALLER)
287 out << " L";
288 else if (*r == CGAL::EQUAL)
289 out << " E";
290 else if (*r == CGAL::LARGER)
291 out << " G";
292 else
293 CGAL_qpe_assertion_msg(false, "incorrect row-type");
294 out << " c" << i << "\n"; // row name is CI
295 }
296
297 // COLUMNS section:
298 typename P::A_iterator a = p.get_a();
299 typename P::C_iterator c = p.get_c();
300 typedef
301 typename std::iterator_traits<typename P::C_iterator>::value_type IT;
302 out << "COLUMNS\n";
303 for (int j=0; j<n; ++j, ++c, ++a) {
304 // make sure that variable appears here even if it has only
305 // zero coefficients
306 bool written = false;
307 if (!CGAL_NTS is_zero (*c)) {
308 out << " x" << j << " obj " << *c << "\n";
309 written = true;
310 }
311 for (int i=0; i<m; ++i) {
312 if (!CGAL_NTS is_zero (*((*a)+i))) {
313 out << " x" << j << " c" << i << " " << *((*a)+i) << "\n";
314 written = true;
315 }
316 }
317 if (!written)
318 out << " x" << j << " obj " << IT(0) << "\n";
319 }
320
321 // RHS section:
322 typename P::B_iterator b = p.get_b();
323 out << "RHS\n";
324 if (!CGAL_NTS is_zero (p.get_c0()))
325 out << " rhs obj " << -p.get_c0() << "\n";
326 for (int i=0; i<m; ++i, ++b)
327 if (!CGAL_NTS is_zero (*b))
328 out << " rhs c" << i << " " << *b << "\n";
329
330 // BOUNDS section:
331 QP_functions_detail::print_bounds (out, p, is_nonnegative);
332
333 // QMATRIX section:
334 QP_functions_detail::print_qmatrix (out, p, is_linear);
335
336 // output end:
337 out << "ENDATA\n";
338 }
339
340 template <typename Program, typename ET,
341 typename Is_linear,typename Is_nonnegative >
solve_program(const Program & p,const ET &,Is_linear,Is_nonnegative,const Quadratic_program_options & options)342 Quadratic_program_solution<ET> solve_program
343 (const Program &p, const ET&,
344 Is_linear,
345 Is_nonnegative,
346 const Quadratic_program_options& options)
347 {
348 typedef QP_solver<
349 Program, ET,
350 QP_solver_impl::QP_tags<Is_linear, Is_nonnegative> >
351 Solver;
352 const Solver* s = new Solver(p, options);
353 Quadratic_program_solution<ET> solution(s);
354 if (options.get_auto_validation()) {
355 // validate solution
356 if (options.get_verbosity() > 0)
357 std::cout << "Validating solution...\n";
358 if (!solution.solves_program(p, Is_linear(), Is_nonnegative())) {
359 // write log file
360 std::ofstream out("QP_solver.log");
361 out << "Error: Program solution is invalid\n"
362 << " (error message: " << solution.get_error() << ")\n"
363 << "------------------\n"
364 << "Solution function:\n"
365 << "------------------\n";
366 print_solution_function (out, Is_linear(), Is_nonnegative());
367 out << "\n"
368 << "--------\n"
369 << "Program:\n"
370 << "--------\n";
371 print_program (out, p, "unsolved", Is_linear(), Is_nonnegative());
372 out << "--------\n"
373 << "Options:\n"
374 << "--------\n"
375 << options << std::endl;
376 // print warning
377 std::cerr
378 << "Error: Program solution is invalid "
379 << "(see QP_solver.log for details)" << std::endl;
380 }
381 }
382 return solution;
383
384 }
385 }
386
387 template <typename QuadraticProgram, typename ET>
solve_quadratic_program(const QuadraticProgram & qp,const ET & dummy,const Quadratic_program_options & options)388 Quadratic_program_solution<ET> solve_quadratic_program
389 (const QuadraticProgram &qp, const ET& dummy,
390 const Quadratic_program_options& options)
391 {
392 return QP_functions_detail::
393 solve_program(qp, dummy, Tag_false(), Tag_false(), options);
394 }
395
396 template <typename NonnegativeQuadraticProgram, typename ET>
solve_nonnegative_quadratic_program(const NonnegativeQuadraticProgram & qp,const ET & dummy,const Quadratic_program_options & options)397 Quadratic_program_solution<ET> solve_nonnegative_quadratic_program
398 (const NonnegativeQuadraticProgram &qp, const ET& dummy,
399 const Quadratic_program_options& options)
400 {
401 return QP_functions_detail::
402 solve_program(qp, dummy, Tag_false(), Tag_true(), options);
403 }
404
405 template <typename LinearProgram, typename ET>
solve_linear_program(const LinearProgram & lp,const ET & dummy,const Quadratic_program_options & options)406 Quadratic_program_solution<ET> solve_linear_program
407 (const LinearProgram &lp, const ET& dummy,
408 const Quadratic_program_options& options)
409 {
410 return QP_functions_detail::
411 solve_program(lp, dummy, Tag_true(), Tag_false(), options);
412 }
413
414 template <typename NonnegativeLinearProgram, typename ET>
solve_nonnegative_linear_program(const NonnegativeLinearProgram & lp,const ET & dummy,const Quadratic_program_options & options)415 Quadratic_program_solution<ET> solve_nonnegative_linear_program
416 (const NonnegativeLinearProgram &lp, const ET& dummy,
417 const Quadratic_program_options& options)
418 {
419 return QP_functions_detail::
420 solve_program(lp, dummy, Tag_true(), Tag_true(), options);
421 }
422
423 } //namespace CGAL
424
425 #endif // CGAL_QP_FUNCTIONS_IMPL_H
426