1 /*-------------------------------------------------------------------------------------*/
2 /* NOMAD - Nonlinear Optimization by Mesh Adaptive Direct search - version 3.7.2 */
3 /* */
4 /* Copyright (C) 2001-2015 Mark Abramson - the Boeing Company, Seattle */
5 /* Charles Audet - Ecole Polytechnique, Montreal */
6 /* Gilles Couture - Ecole Polytechnique, Montreal */
7 /* John Dennis - Rice University, Houston */
8 /* Sebastien Le Digabel - Ecole Polytechnique, Montreal */
9 /* Christophe Tribes - Ecole Polytechnique, Montreal */
10 /* */
11 /* funded in part by AFOSR and Exxon Mobil */
12 /* */
13 /* Author: Sebastien Le Digabel */
14 /* */
15 /* Contact information: */
16 /* Ecole Polytechnique de Montreal - GERAD */
17 /* C.P. 6079, Succ. Centre-ville, Montreal (Quebec) H3C 3A7 Canada */
18 /* e-mail: nomad@gerad.ca */
19 /* phone : 1-514-340-6053 #6928 */
20 /* fax : 1-514-340-5665 */
21 /* */
22 /* This program is free software: you can redistribute it and/or modify it under the */
23 /* terms of the GNU Lesser General Public License as published by the Free Software */
24 /* Foundation, either version 3 of the License, or (at your option) any later */
25 /* version. */
26 /* */
27 /* This program is distributed in the hope that it will be useful, but WITHOUT ANY */
28 /* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A */
29 /* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. */
30 /* */
31 /* You should have received a copy of the GNU Lesser General Public License along */
32 /* with this program. If not, see <http://www.gnu.org/licenses/>. */
33 /* */
34 /* You can find information on the NOMAD software at www.gerad.ca/nomad */
35 /*-------------------------------------------------------------------------------------*/
36 /**
37 \file Pareto_Front.cpp
38 \brief Pareto front (implementation)
39 \author Sebastien Le Digabel
40 \date 2010-04-09
41 \see Pareto_Front.hpp
42 */
43 #include "Pareto_Front.hpp"
44
45 /*------------------------------------------------------*/
46 /* insertion of a point */
47 /* (returns true if the point is a new Pareto point) */
48 /*------------------------------------------------------*/
insert(const NOMAD::Eval_Point & x)49 bool NOMAD::Pareto_Front::insert ( const NOMAD::Eval_Point & x )
50 {
51
52 NOMAD::Pareto_Point pp ( &x );
53
54 if ( _pareto_pts.empty() ) {
55 _pareto_pts.insert (pp);
56 return true;
57 }
58
59 bool insert = false;
60
61 std::set<NOMAD::Pareto_Point>::iterator it = _pareto_pts.begin();
62 while ( it != _pareto_pts.end() ) {
63 if ( pp.dominates (*it) ) {
64 _pareto_pts.erase(it++);
65 insert = true;
66
67 continue;
68 }
69 ++it;
70 }
71
72 if ( !insert ) {
73 insert = true;
74 std::set<NOMAD::Pareto_Point>::iterator end = _pareto_pts.end();
75 for ( it = _pareto_pts.begin() ; it != end ; ++it ) {
76 if ( it->dominates (pp) ) {
77 insert = false;
78 break;
79 }
80 }
81 }
82
83 if ( insert ) {
84 _pareto_pts.insert ( pp );
85 return true;
86 }
87 return false;
88 }
89
90 /*-------------------------------------------------------*/
91 /* get the point that minimizes f2(x) */
92 /* (this is simply the last point of the Pareto front) */
93 /*-------------------------------------------------------*/
get_best_f2(void) const94 const NOMAD::Eval_Point * NOMAD::Pareto_Front::get_best_f2 ( void ) const
95 {
96 if ( _pareto_pts.empty() )
97 return NULL;
98
99 std::set<NOMAD::Pareto_Point>::const_iterator it = _pareto_pts.end();
100 --it;
101
102 return it->get_element();
103 }
104
105 /*------------------------------------------------------*/
106 /* get the reference point */
107 /*------------------------------------------------------*/
get_ref(const NOMAD::Pareto_Point * & xj,NOMAD::Double & delta_j) const108 NOMAD::Point * NOMAD::Pareto_Front::get_ref ( const NOMAD::Pareto_Point *& xj ,
109 NOMAD::Double & delta_j ) const
110 {
111 xj = NULL;
112 delta_j.clear();
113
114 int p = static_cast<int>(_pareto_pts.size());
115
116 // no points in the front:
117 if ( p == 0 )
118 return NULL;
119
120 // just one point in the front:
121 if ( p == 1 ) {
122 xj = &(*_pareto_pts.begin());
123 delta_j = 1.0 / ( xj->get_w() + 1 ); // delta=1.0
124 return NULL;
125 }
126
127 std::set<NOMAD::Pareto_Point>::const_iterator it = _pareto_pts.begin();
128 NOMAD::Point * ref = new NOMAD::Point ( 2 );
129
130 NOMAD::Double f1xm1; // f_1 ( x_{j-1} )
131 NOMAD::Double f1x; // f_1 ( x_j )
132 NOMAD::Double f1xp1; // f_1 ( x_{j+1} )
133
134 NOMAD::Double f2xm1; // f_2 ( x_{j-1} )
135 NOMAD::Double f2x; // f_2 ( x_j )
136 NOMAD::Double f2xp1; // f_2 ( x_{j+1} )
137
138 // two points in the front:
139 if ( p == 2 ) {
140
141 f1xm1 = it->get_f1();
142 f2xm1 = it->get_f2();
143
144 ++it;
145 xj = &(*it);
146
147 f1x = xj->get_f1();
148 f2x = xj->get_f2();
149
150 delta_j = ( (f1x-f1xm1).pow2() + (f2x-f2xm1).pow2() ) / ( xj->get_w() + 1.0 );
151
152 const_cast<NOMAD::Pareto_Point *>(xj)->update_w();
153
154 (*ref)[0] = f1x;
155 (*ref)[1] = f2xm1;
156
157 return ref;
158 }
159
160 // more than two points in the front:
161 std::set<NOMAD::Pareto_Point>::const_iterator end = _pareto_pts.end();
162
163 const NOMAD::Pareto_Point * prev , * cur , * next;
164
165 NOMAD::Double delta;
166
167 prev = &(*it);
168 ++it;
169
170 while ( true ) {
171
172 cur = &(*it);
173
174 ++it;
175 if ( it == end )
176 break;
177
178 next = &(*it);
179
180 f1xm1 = prev->get_f1();
181 f2xm1 = prev->get_f2();
182
183 f1x = cur->get_f1();
184 f2x = cur->get_f2();
185
186 f1xp1 = next->get_f1();
187 f2xp1 = next->get_f2();
188
189 delta = ( (f1x-f1xm1).pow2() + (f2x-f2xm1).pow2() +
190 (f1x-f1xp1).pow2() + (f2x-f2xp1).pow2() ) / ( cur->get_w() + 1.0 );
191
192 if ( !delta_j.is_defined() || delta > delta_j ) {
193 xj = cur;
194 delta_j = delta;
195 (*ref)[0] = f1xp1;
196 (*ref)[1] = f2xm1;
197 }
198
199 prev = cur;
200 }
201
202 const_cast<Pareto_Point *>(xj)->update_w();
203
204 return ref;
205 }
206
207 /*------------------------------------------------------*/
208 /* compute delta and surf */
209 /*------------------------------------------------------*/
get_delta_surf(NOMAD::Double & delta_j,NOMAD::Double & surf,const NOMAD::Point & f_bounds) const210 void NOMAD::Pareto_Front::get_delta_surf ( NOMAD::Double & delta_j ,
211 NOMAD::Double & surf ,
212 const NOMAD::Point & f_bounds ) const
213 {
214 bool def = f_bounds.is_complete();
215 NOMAD::Double f1_min , f1_max , f2_min , f2_max;
216
217 if ( def ) {
218
219 if ( f_bounds.size() == 4 ) {
220
221 f1_min = f_bounds[0];
222 f1_max = f_bounds[1];
223 f2_min = f_bounds[2];
224 f2_max = f_bounds[3];
225
226 if ( f1_min >= f1_max || f2_min >= f2_max ) {
227 f1_min.clear();
228 f1_max.clear();
229 f2_min.clear();
230 f2_max.clear();
231 def = false;
232 }
233 }
234 else
235 def = false;
236 }
237
238 delta_j.clear();
239 surf.clear();
240
241 int p = static_cast<int> ( _pareto_pts.size() );
242
243 // no point in the front:
244 if ( p == 0 ) {
245 if ( def )
246 surf = 1.0;
247 return;
248 }
249
250 const NOMAD::Pareto_Point * xj;
251 NOMAD::Double f1x; // f_1 ( x_j )
252 NOMAD::Double f2x; // f_2 ( x_j )
253
254 NOMAD::Double surf_frame = (def) ?
255 ( f2_max - f2_min ) * ( f1_max - f1_min )
256 : NOMAD::Double();
257
258 // just one point in the front:
259 if ( p == 1 ) {
260 xj = &(*_pareto_pts.begin());
261 delta_j = 1.0 / ( xj->get_w() + 1 ); // delta=1.0
262 f1x = xj->get_f1();
263 f2x = xj->get_f2();
264
265 if ( !def || f1x > f1_max || f1x < f1_min || f2x > f2_max || f2x < f2_min )
266 return;
267
268 surf = ( ( f2_max - f2_min ) * ( f1x - f1_min )
269 + ( f2x - f2_min ) * ( f1_max - f1x ) ) / surf_frame;
270
271 return;
272 }
273
274 std::set<NOMAD::Pareto_Point>::const_iterator it = _pareto_pts.begin();
275
276 NOMAD::Double f1xm1; // f_1 ( x_{j-1} )
277 NOMAD::Double f1xp1; // f_1 ( x_{j+1} )
278
279 NOMAD::Double f2xm1; // f_2 ( x_{j-1} )
280 NOMAD::Double f2xp1; // f_2 ( x_{j+1} )
281
282 // two points in the front:
283 if ( p == 2 ) {
284
285 f1xm1 = it->get_f1();
286 f2xm1 = it->get_f2();
287
288 if ( def && ( f1xm1 < f1_min ||
289 f1xm1 > f1_max ||
290 f2xm1 < f2_min ||
291 f2xm1 > f2_max ) )
292 def = false;
293
294 ++it;
295 xj = &(*it);
296
297 f1x = xj->get_f1();
298 f2x = xj->get_f2();
299
300 if ( def && ( f1x < f1_min ||
301 f1x > f1_max ||
302 f2x < f2_min ||
303 f2x > f2_max ) )
304 def = false;
305
306 delta_j = ( (f1x-f1xm1).pow2() + (f2x-f2xm1).pow2() ) / ( xj->get_w() + 1.0 );
307
308 if ( def )
309 surf = ( ( f2xm1 - f2_min ) * ( f1x - f1xm1 )
310 + ( f2_max - f2_min ) * ( f1xm1 - f1_min )
311 + ( f2x - f2_min ) * ( f1_max - f1x ) ) / surf_frame;
312
313 return;
314 }
315
316 // more than two points in the front:
317 std::set<NOMAD::Pareto_Point>::const_iterator end = _pareto_pts.end();
318
319 const NOMAD::Pareto_Point * prev , * cur , * next;
320
321 NOMAD::Double delta;
322
323 prev = &(*it);
324 f1xm1 = prev->get_f1();
325 f2xm1 = prev->get_f2();
326
327 ++it;
328
329 cur = &(*it);
330 f1x = cur->get_f1();
331
332 if ( def && ( f1xm1 < f1_min ||
333 f1xm1 > f1_max ||
334 f2xm1 < f2_min ||
335 f2xm1 > f2_max ||
336 f1x < f1_min ||
337 f1x > f1_max ) )
338 def = false;
339
340 if ( def )
341 surf = ( f2xm1 - f2_min ) * ( f1x - f1xm1 )
342 + ( f2_max - f2_min ) * ( f1xm1 - f1_min );
343
344 while ( true ) {
345
346 cur = &(*it);
347
348 ++it;
349 if ( it == end )
350 break;
351
352 next = &(*it);
353
354 f1xm1 = prev->get_f1();
355 f2xm1 = prev->get_f2();
356
357 f1x = cur->get_f1();
358 f2x = cur->get_f2();
359
360 f1xp1 = next->get_f1();
361 f2xp1 = next->get_f2();
362
363
364 if ( def &&
365 ( f1xm1 < f1_min || f1xm1 > f1_max || f2xm1 < f2_min || f2xm1 > f2_max ||
366 f1x < f1_min || f1x > f1_max || f2x < f2_min || f2x > f2_max ||
367 f1xp1 < f1_min || f1xp1 > f1_max || f2xp1 < f2_min || f2xp1 > f2_max ) )
368 def = false;
369
370 delta = ( (f1x-f1xm1).pow2() + (f2x-f2xm1).pow2() +
371 (f1x-f1xp1).pow2() + (f2x-f2xp1).pow2() ) / ( cur->get_w() + 1.0 );
372
373 if ( !delta_j.is_defined() || delta > delta_j )
374 delta_j = delta;
375
376 if ( def )
377 surf += ( f2x - f2_min ) * ( f1xp1 - f1x );
378
379 prev = cur;
380 }
381
382 if ( def ) {
383 surf += ( f2xp1 - f2_min ) * ( f1_max - f1xp1 );
384 surf = surf / surf_frame;
385 }
386 else
387 surf.clear();
388 }
389
390 /*------------------------------------------------------*/
391 /* display the Pareto points */
392 /*------------------------------------------------------*/
display(const NOMAD::Display & out) const393 void NOMAD::Pareto_Front::display ( const NOMAD::Display & out ) const
394 {
395 size_t nb = _pareto_pts.size();
396 int cnt = 0;
397 std::set<NOMAD::Pareto_Point>::const_iterator it , end = _pareto_pts.end();
398 for ( it = _pareto_pts.begin() ; it != end ; ++it ) {
399 out << "#";
400 out.display_int_w ( cnt++ , static_cast<int>(nb) );
401 out << " ";
402 it->display ( out );
403 out << std::endl;
404 }
405 }
406
407 /*------------------------------------------------------------------*/
408 /* . begin() and next() methods, to browse the front */
409 /* */
410 /* . example: */
411 /* */
412 /* const Eval_Point * cur = pareto_front.begin(); */
413 /* while (cur) { */
414 /* ... */
415 /* cur = pareto_front.next(); */
416 /* } */
417 /*------------------------------------------------------------------*/
418
419 // begin():
420 // --------
begin(void) const421 const NOMAD::Eval_Point * NOMAD::Pareto_Front::begin ( void ) const
422 {
423 if ( _pareto_pts.empty() )
424 return NULL;
425 _it = _pareto_pts.begin();
426 return _it->get_element();
427 }
428
429 // next(): (supposes that begin() has been called)
430 // -------
next(void) const431 const NOMAD::Eval_Point * NOMAD::Pareto_Front::next ( void ) const
432 {
433 if ( _pareto_pts.empty() )
434 return NULL;
435 ++_it;
436 if ( _it == _pareto_pts.end() )
437 return NULL;
438 return _it->get_element();
439 }
440