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__filtered_base_impl.h $
7 // $Id: QP__filtered_base_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) : Sven Schoenherr
12 // Bernd Gaertner <gaertner@inf.ethz.ch>
13 // Franz Wessendorp
14 // Kaspar Fischer
15
16 namespace CGAL {
17
18 // =============================
19 // class implementation (cont'd)
20 // =============================
21
22 // set-up
23 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
24 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
set()25 set( )
26 {
27 // reserve memory for NT versions of current solution
28 //int l = (std::min)( this->solver().number_of_variables(),
29 // this->solver().number_of_constraints());
30 int l = this->solver().get_l();
31 lambda_NT.resize( l, nt0);
32 set( l, Is_linear());
33 }
34
35 // initialization; BG: who calls this???
36 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
37 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
init()38 init( )
39 {
40 // get properties of quadratic program
41 n = this->solver().number_of_variables();
42 int m = this->solver().number_of_constraints();
43 int s = this->solver().number_of_slack_variables();
44
45 // clear row and column maxima, if necessary
46 if ( ! row_max_A.empty()) row_max_A.clear();
47 if ( ! handled_A.empty()) handled_A.clear();
48 if ( ! col_max .empty()) col_max .clear();
49 init( Is_linear());
50
51 // initialize row and column maxima
52 // assuming nonneg coefficients
53 row_max_c = nt0;
54 C_auxiliary_iterator v_it;
55 for ( v_it = this->solver().c_auxiliary_value_iterator_begin();
56 v_it != this->solver().c_auxiliary_value_iterator_end(); ++v_it) {
57 if (*v_it > row_max_c) row_max_c = (*v_it);
58 }
59
60 handled_A.insert( handled_A.end(), m, false);
61 row_max_A.insert( row_max_A.end(), m, nt0);
62
63 col_max.insert( col_max.end(), n, nt0); // original
64 col_max.insert( col_max.end(), s, nt1); // slack
65 //auxiliary
66 for ( v_it = this->solver().c_auxiliary_value_iterator_begin();
67 v_it != this->solver().c_auxiliary_value_iterator_end(); ++v_it) {
68 col_max.insert( col_max.end(), (*v_it));
69 }
70 }
71
72 // operations
73 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
74 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
init_NT()75 init_NT( )
76 {
77 // ToDo: scale 'x_B_O', 'lambda', and 'd' if necessary
78
79 // get inexact version of 'lambda'
80 std::transform( this->solver().get_lambda_begin(),
81 this->solver().get_lambda_end(),
82 lambda_NT.begin(), et2nt_obj);
83
84 // get inexact version of 'x_B_O'
85 init_NT( Is_linear());
86
87 // get inexact version of 'd'
88 d_NT = et2nt_obj( this->solver().variables_common_denominator());
89 }
90
91 // Here is how it works (Sven's thesis, page 99):
92 // ----------------------------------------------
93 // R_0 := max_j |c_j| (largest objective function coefficient)
94 // R_i^A := max_j |A_ij| (largest entry in row i of A)
95 // R_i^D := max_j 2|D_ij| (largest entry in row i of D)
96 // C_j := max_j (|c_j|, max_i |A_ij|, max_i |D_ij|)
97 //
98 // U := max(
99 // d_NT * R_0,
100 // max_{i basic} |lambda_NT[i]| * R_i^A,
101 // max_{j basic} |x_NT[j]| * R_j^D
102 // )
103 // W := max(
104 // d_NT,
105 // max_{i basic} |lambda_NT[i]|,
106 // max_{j basic} |x_NT[j]|
107 // )
108 // Note: all max_{j basic} are over basic ORIGINAL variables only
109 //
110 // q := (1+1/64) *
111 // (#basic constraints + #basic original variables + 1) *
112 // (#basic constraints + #basic original variables + 2) * 2^(-53)
113 //
114 // the actual error bound for mu_j_NT is then
115 //
116 // min (U * q, W * q * C_j)
117 //
118 // the code maintains and manipulates:
119 //
120 // row_max_c = R_0 (set in init() for phase I, transition() for phase II)
121 // row_max_A[i] = R_i^A (set in update_maxima())
122 // row_max_D[i] = R_i^D (set in update_maxima(Tag_false))
123 // col_max[j] = C_j (set in update_maxima() + Tag_false variant)
124 // bound1 = U * q (computed in update_maxima() + Tag_true/Tag_false )
125 // bound2_wq = W * q (computed in update_maxima() + Tag_true/Tag_false )
126 //
127 // Changes for upper-bounded case:
128 // -------------------------------
129 // in the scalar product considered on p.100, we have another term:
130 // d * w_j, where w_j is the correction term 2 x_N^T D_Nj that shows
131 // up additionally in the pricing and that the QP_solver maintains in
132 // the vector w. The first bound on p.101 yields U, while the second
133 // bound yields W. This means that
134 // - C_j = max_y|y_i| has an additional term |w_NT[j]| in the maximum,
135 // - U has an additional term d_NT * R_w in the maximum, where
136 //
137 // R_w = max_j |w_j|
138 //
139 // The problem is that R_w is not a static quantity like max_j|c_j|,
140 // so it has to be handled per j in certify_mu_j_NT:
141 // --> for q: (c+b)*(c+b+1) -> (c+b+1)*(c+b+2)
142 // --> for bound1: take d_NT * |w_NT[j]| into account for U
143 // --> for bound2: take |w_NT[j]| into account for C_j
144
145
146 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
147 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
update_maxima()148 update_maxima( )
149 {
150 // get properties of quadratic program
151 int c = this->solver().number_of_basic_constraints();
152 int b = this->solver().number_of_basic_original_variables();
153
154 // initialize error bounds
155 bound1 = d_NT * row_max_c;
156 bound2_wq = d_NT;
157
158 // update row and column maxima of 'A'
159 A_iterator a_it = this->solver().a_begin();
160 R_iterator r_it = this->solver().row_type_begin();
161 int row, col;
162 NT row_max, z;
163
164 Basic_constraint_index_iterator it;
165 Values_NT_iterator v_it = lambda_NT.begin();
166 for ( it = this->solver().basic_constraint_indices_begin();
167 it != this->solver().basic_constraint_indices_end(); ++it, ++v_it) {
168 row = *it;
169
170 // row not handled yet?
171 if ( ! handled_A[ row]) {
172
173 // slack variable (or artificial in phase I) involved?
174 row_max = ( ( r_it[ row] != CGAL::EQUAL)
175 || ( this->solver().phase() == 1) ? nt1 : nt0);
176
177 // scan row and update maxima
178 for ( col = 0; col < n; ++col) {
179 z = CGAL::abs( *((*(a_it + col))+row));
180 if ( z > row_max ) row_max = z;
181 if ( z > col_max[ col]) col_max[ col] = z;
182 }
183 row_max_A[ row] = row_max;
184 handled_A[ row] = true;
185 }
186 // update bounds
187 z = CGAL::abs( *v_it);
188 if ( z > bound2_wq) bound2_wq = z;
189 z *= row_max_A[ row];
190 if ( z > bound1) bound1 = z;
191 }
192
193 // update row and column maxima of 'D', if necessary
194 if (this->solver().phase() == 1) {
195 update_maxima( Tag_true());
196 } else {
197 update_maxima( Is_linear());
198 }
199
200 // finalize error bounds
201 // ToDo: use std::numeric_limits for 'machine epsilon' to
202 // support types different from double
203 set_q(c, b);
204 bound1 *= q;
205 bound2_wq *= q;
206
207 CGAL_qpe_debug {
208 this->vout() << std::endl
209 << "first bound for certification: " << bound1 << std::endl;
210 }
211 }
212
213 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ > // QP case
214 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
update_maxima(Tag_false)215 update_maxima( Tag_false)
216 {
217 // update row and column maxima of 'D'
218 D_row_iterator d_row_it;
219 int row, col;
220 NT row_max, z;
221
222 Basic_variable_index_iterator it;
223 Values_NT_iterator v_it = x_B_O_NT.begin();
224 for ( it = this->solver().basic_original_variable_indices_begin();
225 it != this->solver().basic_original_variable_indices_end();
226 ++it, ++v_it) {
227 row = *it;
228
229 // row not handled yet?
230 if ( ! handled_D[ row]) {
231 d_row_it = this->solver().d_begin()[ row];
232
233 // scan row and update maxima
234 row_max = nt0;
235 for ( col = 0; col < n; ++col, ++d_row_it) {
236 z = CGAL::abs( *d_row_it);
237 if ( z > row_max ) row_max = z;
238 if ( z > col_max[ col]) col_max[ col] = z;
239 }
240 row_max_D[ row] = row_max;
241 handled_D[ row] = true;
242 }
243 // update bounds
244 z = CGAL::abs( (*v_it));
245 if ( z > bound2_wq) bound2_wq = z;
246 z *= row_max_D[ row];
247 if ( z > bound1) bound1 = z;
248 }
249 }
250
251 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ > // LP case
252 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
update_maxima(Tag_true)253 update_maxima( Tag_true)
254 {
255
256 // update basic original variables maxima
257 NT z;
258
259 Values_NT_iterator v_it;
260 for ( v_it = x_B_O_NT.begin(); v_it != x_B_O_NT.end(); ++v_it) {
261
262 // update bounds
263 z = CGAL::abs( (*v_it));
264 if ( z > bound2_wq) bound2_wq = z;
265 }
266 }
267
268 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
269 bool QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
certify_mu_j_NT(int j,Tag_true)270 certify_mu_j_NT( int j, Tag_true) const // standard form
271 {
272 // compute 'mu_j' with inexact arithmetic again
273 NT mu = mu_j_NT( j);
274
275 CGAL_qpe_debug {
276 this->vout() << "mu_" << j << " [NT]: " << mu;
277 }
278
279 // check against first bound
280 if ( mu >= bound1) {
281 CGAL_qpe_debug {
282 this->vout() << " ok [1]" << std::endl;
283 }
284 return true;
285 }
286
287 // compute and check against second bound
288 NT bound2 = bound2_wq * col_max[ j];
289 if ( mu >= bound2) {
290 CGAL_qpe_debug {
291 this->vout() << " ok [2: " << bound2 << "]" << std::endl;
292 }
293 return true;
294 }
295
296 // compute and check exact 'mu_j'
297 ET mu_et = this->mu_j( j);
298
299 CGAL_qpe_debug {
300 this->vout() << " " << ( mu_et >= this->et0 ? "ok" : "MISSED")
301 << " [exact: " << mu_et << "]" << std::endl;
302 }
303 return ( mu_et >= this->et0);
304 }
305
306 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
307 bool QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
certify_mu_j_NT(int j,Tag_false)308 certify_mu_j_NT( int j, Tag_false) const // case with bounds
309 {
310 // compute 'mu_j' with inexact arithmetic again; this call sets w_j_NT
311 NT mu = CGAL::abs(mu_j_NT( j)); // may have both signs in bounded case
312 NT w_j_NT_abs = CGAL::abs(w_j_NT);
313
314 CGAL_qpe_debug {
315 this->vout() << "|mu_|" << j << " [NT]: " << mu;
316 }
317
318 // check against first bound, taking w_j_NT into account
319 // q has been set in previous call to update_maxima()
320 // bound1_with_w = (U + d_NT*|w_j|) * q
321 NT bound1_with_w = bound1 + (d_NT * w_j_NT_abs * q);
322 if ( mu >= bound1_with_w) {
323 CGAL_qpe_debug {
324 this->vout() << " ok [1]" << std::endl;
325 }
326 return true;
327 }
328
329 // compute and check against second bound, taking w_j_NT into account
330 // bound2_with_w = max (W, |w_j|)
331 NT bound2_with_w = bound2_wq * (std::max)(col_max[ j], w_j_NT_abs);
332 if ( mu >= bound2_with_w) {
333 CGAL_qpe_debug {
334 this->vout() << " ok [2: " << bound2_with_w << "]" << std::endl;
335 }
336 return true;
337 }
338
339 // the bounds are insufficient: compute and check exact 'mu_j'
340 ET mu_et = this->mu_j( j);
341 bool improving = is_improving(j, mu_et, this->et0);
342
343 CGAL_qpe_debug {
344 this->vout() << " " << (!improving ? "ok" : "MISSED")
345 << " [exact: " << mu_et << "]" << std::endl;
346 }
347 return ( !improving);
348 }
349
350 // transition
351 template < typename Q, typename ET, typename Tags, class NT_, class ET2NT_ >
352 void QP__filtered_base<Q,ET,Tags,NT_,ET2NT_>::
transition()353 transition( )
354 {
355
356 // update row and column maxima with original objective vector 'c'
357 C_iterator c_it = this->solver().c_begin();
358 NT z;
359 row_max_c = nt0;
360 for ( int i = 0; i < n; ++i, ++c_it) {
361 z = CGAL::abs( *c_it);
362 if ( z > row_max_c ) row_max_c = z;
363 if ( z > col_max[ i]) col_max[ i] = z;
364 }
365
366 // initialize row maxima of 'D', if necessary
367 transition( n, Is_linear());
368 }
369
370 } //namespace CGAL
371
372 // ===== EOF ==================================================================
373