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 TGP_Model_Search.cpp
38 \brief TGP Model search (implementation)
39 \author Sebastien Le Digabel
40 \date 2011-02-17
41 \see TGP_Model_Search.hpp
42 */
43
44 #ifndef USE_TGP
45
46 int TGP_MODEL_SEARCH_DUMMY; // avoids that TGP_Model_Search.o has no symbols with ranlib
47
48 #else
49
50 #include "TGP_Model_Search.hpp"
51
52 /*-----------------------------------*/
53 /* reset (virtual) */
54 /*-----------------------------------*/
reset(void)55 void NOMAD::TGP_Model_Search::reset ( void )
56 {
57 if ( _model )
58 delete _model;
59 _model = NULL;
60 }
61
62 /*--------------------------------------------------------*/
63 /* delete a list of points: one version for points, and */
64 /* one version for evaluation points (static, private) */
65 /*--------------------------------------------------------*/
clear_pts(std::vector<NOMAD::Point * > & pts)66 void NOMAD::TGP_Model_Search::clear_pts ( std::vector<NOMAD::Point *> & pts )
67 {
68 size_t k , n = pts.size();
69 for ( k = 0 ; k < n ; ++k )
70 delete pts[k];
71 pts.clear();
72 }
73
clear_pts(std::vector<NOMAD::Eval_Point * > & pts)74 void NOMAD::TGP_Model_Search::clear_pts ( std::vector<NOMAD::Eval_Point *> & pts )
75 {
76 size_t k , n = pts.size();
77 for ( k = 0 ; k < n ; ++k )
78 delete pts[k];
79 pts.clear();
80 }
81
82 /*------------------------------------------------------------------*/
83 /* the search */
84 /*------------------------------------------------------------------*/
85 /* Search parameters: */
86 /* ------------------ */
87 /* */
88 /* . MODEL_SEARCH: flag to activate the model search (MS) */
89 /* (here its value is NOMAD::TGP_MODEL) */
90 /* */
91 /* . MODEL_SEARCH_OPTIMISTIC: if true, the direction from the */
92 /* model center to the trial point */
93 /* is computed and prossibly used */
94 /* in the speculative search */
95 /* default=yes */
96 /* */
97 /* . MODEL_SEARCH_PROJ_TO_MESH: project or not to mesh */
98 /* default=yes */
99 /* */
100 /* . MODEL_SEARCH_MAX_TRIAL_PTS: limit on the number of trial */
101 /* points for one search */
102 /* default=10 */
103 /* */
104 /* . MODEL_TGP_MODE: TGP mode (FAST or PRECISE) */
105 /* default=FAST */
106 /* */
107 /*------------------------------------------------------------------*/
search(NOMAD::Mads & mads,int & nb_search_pts,bool & stop,NOMAD::stop_type & stop_reason,NOMAD::success_type & success,bool & count_search,const NOMAD::Eval_Point * & new_feas_inc,const NOMAD::Eval_Point * & new_infeas_inc)108 void NOMAD::TGP_Model_Search::search ( NOMAD::Mads & mads ,
109 int & nb_search_pts ,
110 bool & stop ,
111 NOMAD::stop_type & stop_reason ,
112 NOMAD::success_type & success ,
113 bool & count_search ,
114 const NOMAD::Eval_Point *& new_feas_inc ,
115 const NOMAD::Eval_Point *& new_infeas_inc )
116 {
117 new_feas_inc = new_infeas_inc = NULL;
118 nb_search_pts = 0;
119 success = NOMAD::UNSUCCESSFUL;
120 count_search = false;
121
122 _one_search_stats.reset();
123
124 const NOMAD::Display & out = _p.out();
125 NOMAD::dd_type display_degree = out.get_search_dd();
126 int display_lim = 15;
127
128 if ( stop ) {
129 if ( display_degree == NOMAD::FULL_DISPLAY )
130 out << "TGP_Model_Search::search(): not performed (stop flag is active)"
131 << std::endl;
132 return;
133 }
134
135 // active cache (we accept only true function evaluations):
136 const NOMAD::Cache & cache = mads.get_cache();
137 if ( cache.get_eval_type() != NOMAD::TRUTH ) {
138 if ( display_degree == NOMAD::FULL_DISPLAY )
139 out << "TGP_Model_Search::search(): not performed on surrogates"
140 << std::endl;
141 return;
142 }
143
144 // check that there is one objective exactly:
145 const std::list<int> & index_obj_list = _p.get_index_obj();
146
147 if ( index_obj_list.empty() ) {
148 if ( display_degree == NOMAD::FULL_DISPLAY )
149 out << "TGP_Model_Search::search(): not performed with no objective function"
150 << std::endl;
151 return;
152 }
153 if ( index_obj_list.size() > 1 ) {
154 if ( display_degree == NOMAD::FULL_DISPLAY )
155 out << "TGP_Model_Search::search(): not performed with biobjective"
156 << std::endl;
157 return;
158 }
159
160 // active barrier:
161 const NOMAD::Barrier & barrier = mads.get_true_barrier();
162
163 // get the incumbent:
164 const NOMAD::Eval_Point * incumbent = barrier.get_best_feasible();
165 if ( !incumbent )
166 incumbent = barrier.get_best_infeasible();
167 if ( !incumbent ) {
168 if ( display_degree == NOMAD::FULL_DISPLAY )
169 out << "TGP_Model_Search::search(): no incumbent"
170 << std::endl;
171 return;
172 }
173
174 // get and check the signature, and compute the dimension:
175 NOMAD::Signature * signature = incumbent->get_signature();
176
177 if ( !signature ) {
178 if ( display_degree == NOMAD::FULL_DISPLAY )
179 out << "TGP_Model_Search::search(): no signature"
180 << std::endl;
181 return;
182 }
183
184 int n = signature->get_n();
185
186 if ( n != incumbent->size() ) {
187 if ( display_degree == NOMAD::FULL_DISPLAY )
188 out << "TGP_Model_Search::search(): incompatible signature"
189 << std::endl;
190 return;
191 }
192
193 // initial displays:
194 if ( display_degree == NOMAD::FULL_DISPLAY ) {
195 std::ostringstream oss;
196 oss << NOMAD::MODEL_SEARCH << " #"
197 << _all_searches_stats.get_MS_nb_searches();
198 out << std::endl << NOMAD::open_block ( oss.str() ) << std::endl;
199 }
200
201 // from this point the search is counted:
202 count_search = true;
203 _one_search_stats.add_MS_nb_searches();
204
205 NOMAD::Point delta_m;
206 signature->get_mesh()->get_delta(delta_m);
207
208 // initial displays:
209 if ( display_degree == NOMAD::FULL_DISPLAY ) {
210 #ifdef TGP_DEBUG
211 out << "seed : "
212 << _p.get_seed() << std::endl;
213 #endif
214 out << "number of cache points: " << cache.size() << std::endl
215 << "mesh size parameter : ( " << delta_m << " )" << std::endl
216 << "incumbent : ( ";
217 incumbent->NOMAD::Point::display
218 ( out , " " , 2 , NOMAD::Point::get_display_limit() );
219 out << " )" << std::endl;
220 }
221
222 // construct the model:
223 NOMAD::Stats & stats = mads.get_stats();
224 bool compute_Ds2x;
225 std::vector<NOMAD::Eval_Point *> XX;
226 std::string error_str;
227
228 if ( !model_construction ( cache ,
229 *incumbent ,
230 delta_m ,
231 out ,
232 display_degree ,
233 display_lim ,
234 stats ,
235 compute_Ds2x ,
236 XX ,
237 stop ,
238 stop_reason ,
239 error_str ) ) {
240 if ( display_degree == NOMAD::FULL_DISPLAY )
241 out << NOMAD::close_block ( "failure: " + error_str )
242 << std::endl;
243 return;
244 }
245
246 // trial_pts = oracle_pts + Ds2x_pts + improv_pts
247 //
248 // oracle_pts: given by the model optimization
249 // Ds2x_pts : XX points with the largest expected reduction in predictive variance
250 // improv_pts: XX points with the largest expected improvement for the objective
251
252 int max_pts = _p.get_model_search_max_trial_pts();
253
254 /*-----------------------*/
255 /* oracle points (1/3) */
256 /*-----------------------*/
257 std::vector<NOMAD::Point *> oracle_pts;
258 if ( !create_oracle_pts ( cache ,
259 *incumbent ,
260 delta_m ,
261 out ,
262 display_degree ,
263 display_lim ,
264 XX ,
265 oracle_pts ,
266 stop ,
267 stop_reason ) && stop ) {
268
269 // delete XX and oracle_pts:
270 NOMAD::TGP_Model_Search::clear_pts ( XX );
271 NOMAD::TGP_Model_Search::clear_pts ( oracle_pts );
272
273 // quit:
274 if ( display_degree == NOMAD::FULL_DISPLAY )
275 out << NOMAD::close_block ( "algorithm stop" )
276 << std::endl;
277 return;
278 }
279
280 /*---------------------*/
281 /* Ds2x points (2/3) */
282 /*---------------------*/
283 std::vector<NOMAD::Point *> Ds2x_pts;
284 if ( compute_Ds2x )
285 create_Ds2x_pts ( XX , out , display_degree , display_lim , Ds2x_pts );
286
287 /*-----------------------*/
288 /* improv points (3/3) */
289 /*-----------------------*/
290 std::vector<NOMAD::Point *> improv_pts;
291 create_improv_pts ( XX ,
292 *incumbent ,
293 max_pts ,
294 out ,
295 display_degree ,
296 display_lim ,
297 improv_pts );
298
299 // create the complete list of trial points:
300 // -----------------------------------------
301 std::vector<NOMAD::Point *> trial_pts;
302 create_trial_pts ( oracle_pts ,
303 Ds2x_pts ,
304 improv_pts ,
305 *incumbent ,
306 max_pts ,
307 out ,
308 display_degree ,
309 trial_pts );
310
311 // evaluator control:
312 NOMAD::Evaluator_Control & ev_control = mads.get_evaluator_control();
313
314 // add the trial points to the evaluator control for evaluation:
315 int i , nop = trial_pts.size();
316 for ( i = 0 ; i < nop ; ++i )
317 register_point ( *trial_pts[i] ,
318 *signature ,
319 *incumbent ,
320 display_degree ,
321 ev_control );
322
323 // display the evaluator control list of points:
324 if ( display_degree == NOMAD::FULL_DISPLAY ) {
325 out << std::endl << NOMAD::open_block ( "list of trial points" );
326 const std::set<NOMAD::Priority_Eval_Point> & lop = ev_control.get_eval_lop();
327 std::set<NOMAD::Priority_Eval_Point>::const_iterator it , end = lop.end();
328 nop = lop.size();
329 for ( it = lop.begin() , i = 0 ; it != end ; ++it , ++i ) {
330 out << "#";
331 out.display_int_w ( i , nop );
332 out << " x=( ";
333 it->get_point()->NOMAD::Point::display ( out , " " , 15 , -1 );
334 out << " )" << std::endl;
335 }
336 out.close_block();
337 }
338
339 // delete XX and the trial points
340 // (do not delete Ds2x_pts and improv_pts because
341 // they are XX points, contrary to oracle_pts):
342 NOMAD::TGP_Model_Search::clear_pts ( XX );
343 NOMAD::TGP_Model_Search::clear_pts ( oracle_pts );
344
345 nb_search_pts = ev_control.get_nb_eval_points();
346
347 if ( nb_search_pts == 0 ) {
348 if ( display_degree == NOMAD::FULL_DISPLAY )
349 out << std::endl << "no trial point" << std::endl;
350 }
351
352 else {
353
354 _one_search_stats.update_MS_max_search_pts ( nb_search_pts );
355
356 // evaluate the trial points:
357 // --------------------------
358 int bbe = stats.get_bb_eval();
359 int sgte_eval = stats.get_sgte_eval ();
360 int cache_hits = stats.get_cache_hits();
361
362 new_feas_inc = new_infeas_inc = NULL;
363
364 ev_control.disable_model_eval_sort();
365
366 std::list<const NOMAD::Eval_Point *> * evaluated_pts = NULL;
367 if ( display_degree == NOMAD::FULL_DISPLAY )
368 evaluated_pts = new std::list<const NOMAD::Eval_Point *>;
369
370 ev_control.eval_list_of_points ( _type ,
371 mads.get_true_barrier() ,
372 mads.get_sgte_barrier() ,
373 mads.get_pareto_front() ,
374 stop ,
375 stop_reason ,
376 new_feas_inc ,
377 new_infeas_inc ,
378 success ,
379 evaluated_pts );
380
381 // display the prediction error for the evaluated points:
382 if ( display_degree == NOMAD::FULL_DISPLAY ) {
383 display_eval_pred_errors ( *evaluated_pts , out );
384 delete evaluated_pts;
385 }
386
387 ev_control.enable_model_eval_sort();
388
389 // update stats:
390 _one_search_stats.add_MS_bb_eval ( stats.get_bb_eval () - bbe );
391 _one_search_stats.add_MS_sgte_eval ( stats.get_sgte_eval () - sgte_eval );
392 _one_search_stats.add_MS_cache_hits ( stats.get_cache_hits() - cache_hits );
393
394 if ( success == NOMAD::FULL_SUCCESS )
395 _one_search_stats.add_MS_success();
396
397 _one_search_stats.add_MS_pts ( nb_search_pts );
398 }
399
400 // update stats objects:
401 stats.update_model_stats ( _one_search_stats );
402 _all_searches_stats.update ( _one_search_stats );
403
404 // final display:
405 if ( display_degree == NOMAD::FULL_DISPLAY ) {
406 std::ostringstream oss;
407 oss << "end of " << NOMAD::MODEL_SEARCH << " (" << success << ")";
408 out << std::endl << NOMAD::close_block ( oss.str() ) << std::endl;
409 }
410 }
411
412 /*---------------------------------------------------------------*/
413 /* create XX a list of prediction points (private) */
414 /*---------------------------------------------------------------*/
set_XX(const NOMAD::Cache & cache,int n,int m,const NOMAD::Point & incumbent,const NOMAD::Point & delta_m,std::vector<NOMAD::Eval_Point * > & XX) const415 void NOMAD::TGP_Model_Search::set_XX
416 ( const NOMAD::Cache & cache ,
417 int n ,
418 int m ,
419 const NOMAD::Point & incumbent ,
420 const NOMAD::Point & delta_m ,
421 std::vector<NOMAD::Eval_Point *> & XX ) const
422 {
423 // . we begin with 1999 points, filter them to remove points
424 // that appear more than once or in the cache.
425 // . then we prune the list to 499 points.
426 // . we do not begin directly with 500 points because it
427 // gives more flexibility with the projection.
428 // . we terminate by adding the incumbent.
429 // --> we obtain a list of at most 500 prediction points.
430
431 int j , i = 0 , n_XX = 1999; // 2000-1
432 bool proj_to_mesh = _p.get_model_search_proj_to_mesh();
433 bool remove_pt;
434
435 // XX is made of n_XX-1 LH points inside the model bounds
436 // (and the incumbent will be added at the end):
437 NOMAD::LH_Search::LH_points ( n, m, n_XX, _model->get_lb(), _model->get_ub(), XX );
438
439 while ( i < n_XX ) {
440
441 remove_pt = false;
442
443 // project to the mesh:
444 if ( proj_to_mesh )
445 XX[i]->project_to_mesh ( incumbent , delta_m , _p.get_lb() , _p.get_ub() );
446
447 // remove if point is in cache:
448 if ( cache.find ( *XX[i] ) )
449 remove_pt = true;
450
451 // check if this point is already in XX
452 // (may occur only if the point has been projected):
453 if ( proj_to_mesh && !remove_pt ) {
454
455 if ( incumbent == (*XX[i]) )
456 remove_pt = true;
457
458 else
459 for ( j = 0 ; j < i ; ++j )
460 if ( XX[j]->NOMAD::Point::operator == ( *XX[i] ) ) {
461 remove_pt = true;
462 break;
463 }
464 }
465
466 // remove the point:
467 if ( remove_pt ) {
468 delete XX[i];
469 --n_XX;
470 if ( i != n_XX )
471 XX[i] = XX[n_XX];
472 XX.resize ( n_XX );
473 }
474 else
475 ++i;
476 }
477
478 // reduce to 500-1 points (we eliminate randomly):
479 while ( n_XX >= 500 ) {
480
481 i = rand()%n_XX;
482
483 delete XX[i];
484 --n_XX;
485 if ( i != n_XX )
486 XX[i] = XX[n_XX];
487 XX.resize(n_XX);
488 }
489
490 // add the incumbent as the last point of XX:
491 XX.push_back ( new NOMAD::Eval_Point ( n , m ) );
492 XX[n_XX]->NOMAD::Point::operator = ( incumbent );
493 }
494
495 /*---------------------------------------------------------------------*/
496 /* create the complete list of trial points (oracle + Ds2x + improv) */
497 /* (private) */
498 /*---------------------------------------------------------------------*/
create_trial_pts(const std::vector<NOMAD::Point * > & oracle_pts,const std::vector<NOMAD::Point * > & Ds2x_pts,const std::vector<NOMAD::Point * > & improv_pts,const NOMAD::Point & incumbent,int max_pts,const NOMAD::Display & out,NOMAD::dd_type display_degree,std::vector<NOMAD::Point * > & trial_pts) const499 void NOMAD::TGP_Model_Search::create_trial_pts
500 ( const std::vector<NOMAD::Point *> & oracle_pts , // IN
501 const std::vector<NOMAD::Point *> & Ds2x_pts , // IN
502 const std::vector<NOMAD::Point *> & improv_pts , // IN
503 const NOMAD::Point & incumbent , // IN
504 int max_pts , // IN
505 const NOMAD::Display & out , // IN
506 NOMAD::dd_type display_degree , // IN
507 std::vector<NOMAD::Point *> & trial_pts ) const // OUT
508 {
509 bool found;
510 size_t i , j ,
511 n1 = oracle_pts.size() ,
512 n2 = Ds2x_pts.size () ,
513 n3 = improv_pts.size();
514
515 std::vector<NOMAD::Point *> l2 , l3;
516
517 // 1. remove duplicates:
518 // ---------------------
519 //
520 // . oracle_pts are not XX points
521 // . Ds2x_pts and improv_pts are XX points
522 // . there is no duplicate in each list separately
523 //
524 for ( i = 0 ; i < n2 ; ++i ) {
525 found = false;
526
527 if ( *Ds2x_pts[i] == incumbent )
528 found = true;
529 else {
530 for ( j = 0 ; j < n1 ; ++j )
531 if ( *Ds2x_pts[i] == *oracle_pts[j] ) {
532 found = true;
533 break;
534 }
535 }
536 if ( !found )
537 l2.push_back ( Ds2x_pts[i] );
538 }
539
540 n2 = l2.size();
541
542 for ( i = 0 ; i < n3 ; ++i ) {
543 found = false;
544 if ( *improv_pts[i] == incumbent )
545 found = true;
546 else {
547 for ( j = 0 ; j < n1 ; ++j )
548 if ( *improv_pts[i] == *oracle_pts[j] ) {
549 found = true;
550 break;
551 }
552 }
553 if ( !found ) {
554 for ( j = 0 ; j < n2 ; ++j )
555 if ( improv_pts[i] == l2[j] ) {
556 found = true;
557 break;
558 }
559 }
560 if ( !found )
561 l3.push_back ( improv_pts[i] );
562 }
563
564 n3 = l3.size();
565
566 // 2. construct the list of trial points:
567 // ------------------------------------
568 trial_pts.clear();
569
570 int nb_pts = static_cast<int> ( n1 + n2 + n3 );
571
572 // no need to reduce the number of trial points:
573 if ( max_pts <= 0 || nb_pts <= max_pts ) {
574
575
576 // 1. oracle points:
577 for ( i = 0 ; i < n1 ; ++i )
578 trial_pts.push_back ( oracle_pts[i] );
579
580 // 2, improv points:
581 for ( i = 0 ; i < n3 ; ++i )
582 trial_pts.push_back ( l3[i] );
583
584 // 3. Ds2x points:
585 for ( i = 0 ; i < n2 ; ++i )
586 trial_pts.push_back ( l2[i] );
587
588 }
589
590 // reduce the list to max_pts points:
591 else {
592
593 nb_pts = 0;
594
595 size_t i1 = 0 , i2 = 0 , i3 = 0;
596
597 while ( true ) {
598
599 // one point from the oracle points:
600 if ( i1 < n1 ) {
601 trial_pts.push_back ( oracle_pts[i1++] );
602 ++nb_pts;
603 if ( nb_pts == max_pts )
604 break;
605 }
606
607 // two from the improv points:
608 for ( i = 0 ; i < 2 ; ++i ) {
609 if ( i3 < n3 ) {
610 trial_pts.push_back ( l3[i3++] );
611 ++nb_pts;
612 if ( nb_pts == max_pts )
613 break;
614 }
615 }
616 if ( nb_pts == max_pts )
617 break;
618
619 // one from the Ds2x points:
620 if ( i2 < n2 ) {
621 trial_pts.push_back ( l2[i2++] );
622 ++nb_pts;
623 if ( nb_pts == max_pts )
624 break;
625 }
626 }
627 }
628
629 // 3. display the list of trial points:
630 // ------------------------------------
631 // if ( display_degree == NOMAD::FULL_DISPLAY ) {
632 // out.open_block ( "list of trial points (debug)" );
633 // n1 = trial_pts.size();
634 // for ( i = 0 ; i < n1 ; ++i ) {
635 // out << "#";
636 // out.display_int_w ( i , n1 );
637 // out << " x=( " << *trial_pts[i] << " )" << std::endl;
638 // }
639 // out.close_block();
640 // }
641 }
642
643 /*---------------------------------------------------------*/
644 /* create the list of improv points (private) */
645 /*---------------------------------------------------------*/
create_improv_pts(const std::vector<NOMAD::Eval_Point * > & XX,const NOMAD::Point & incumbent,int max_pts,const NOMAD::Display & out,NOMAD::dd_type display_degree,int display_lim,std::vector<NOMAD::Point * > & improv_pts) const646 void NOMAD::TGP_Model_Search::create_improv_pts
647 ( const std::vector<NOMAD::Eval_Point *> & XX , // IN
648 const NOMAD::Point & incumbent , // IN
649 int max_pts , // IN
650 const NOMAD::Display & out , // IN
651 NOMAD::dd_type display_degree , // IN
652 int display_lim , // IN
653 std::vector<NOMAD::Point *> & improv_pts ) const // OUT
654 {
655 if ( display_degree == NOMAD::FULL_DISPLAY )
656 out << NOMAD::open_block ( "improv points construction" )
657 << std::endl;
658
659 // reset improv points:
660 improv_pts.clear();
661
662 // display improv points directly from the NOMAD::TGP_Model object:
663 #ifdef TGP_DEBUG
664 out << NOMAD::open_block ( "expected improvement of the objective (improv)" );
665 _model->display_improv ( out , display_lim );
666 out << NOMAD::close_block() << std::endl;
667 #endif
668
669 std::list<int> pts_indexes;
670 _model->get_improv_points ( pts_indexes );
671
672 if ( pts_indexes.empty() ) {
673 if ( display_degree == NOMAD::FULL_DISPLAY )
674 out << std::endl
675 << NOMAD::close_block ( "no improv candidate" )
676 << std::endl;
677 return;
678 }
679
680 std::list<int>::const_iterator it , end = pts_indexes.end();
681
682 // with constraints, the list is re-sorted in order to put the
683 // predicred feasible points first:
684 if ( _p.has_constraints() ) {
685
686 std::list<int> feas_pts , infeas_pts;
687
688 NOMAD::Double h , f;
689 const NOMAD::Double & h_min = _p.get_h_min();
690
691 for ( it = pts_indexes.begin() ; it != end ; ++it ) {
692 if ( predict ( *XX[*it] , h , f ) && h <= h_min )
693 feas_pts.push_back ( *it );
694 else
695 infeas_pts.push_back ( *it );
696 }
697
698 pts_indexes.clear();
699
700 end = feas_pts.end();
701 for ( it = feas_pts.begin() ; it != end ; ++it )
702 pts_indexes.push_back ( *it );
703
704 end = infeas_pts.end();
705 for ( it = infeas_pts.begin() ; it != end ; ++it )
706 pts_indexes.push_back ( *it );
707
708 end = pts_indexes.end();
709 }
710
711 // compute max_index just for the display:
712 int i , j , max_index = -1 , ni = -1;
713
714 if ( display_degree == NOMAD::FULL_DISPLAY ) {
715 ni = pts_indexes.size();
716 if ( max_pts > 0 && ni > max_pts )
717 ni = max_pts;
718
719 for ( it = pts_indexes.begin() ; it != end ; ++it ) {
720 if ( *it > max_index )
721 max_index = *it;
722 }
723 }
724
725 // add the points to improv_pts:
726 bool rejected;
727 i = j = 0;
728 for ( it = pts_indexes.begin() ; it != end ; ++it ) {
729
730 // we check the max number of points:
731 rejected = ( max_pts > 0 && max_pts == i );
732
733 // we reject the point if it is the incumbent:
734 if ( !rejected && incumbent == *XX[*it] )
735 rejected = true;
736
737 // we add the point:
738 if ( !rejected ) {
739 improv_pts.push_back ( XX[*it] );
740 ++i;
741 }
742
743 // display:
744 if ( display_degree == NOMAD::FULL_DISPLAY ) {
745 if ( display_lim <= 0 || j < display_lim ) {
746 if ( rejected )
747 out << "rejected candidate ";
748 else {
749 out << "improv candidate #";
750 out.display_int_w ( i-1 , ni );
751 }
752 out << " (XX point #";
753 out.display_int_w ( *it , max_index );
754 out << "): x=( ";
755 XX[*it]->NOMAD::Point::display ( out , " " , 6 , -1 );
756 out << " )";
757 if ( rejected ) {
758 if ( max_pts > 0 && max_pts == i )
759 out << " (max number of points)";
760 else
761 out << " (candidate==incumbent)";
762 }
763 else
764 out << " improv=" << _model->get_improv(*it);
765 out << std::endl;
766 }
767 if ( display_lim > 0 && j == display_lim )
768 out << "..." << std::endl;
769 ++j;
770 }
771
772 // if no display, stop the loop if there is already too many points:
773 else if ( max_pts > 0 && max_pts == i )
774 break;
775 }
776
777 if ( display_degree == NOMAD::FULL_DISPLAY )
778 out << std::endl
779 << NOMAD::close_block ( "end of improv points construction" )
780 << std::endl;
781 }
782
783 /*------------------------------------------------------------*/
784 /* create the list of Ds2x points, the points that maximize */
785 /* the expected reduction in predictive variance */
786 /* (private) */
787 /*------------------------------------------------------------*/
create_Ds2x_pts(const std::vector<NOMAD::Eval_Point * > & XX,const NOMAD::Display & out,NOMAD::dd_type display_degree,int display_lim,std::vector<NOMAD::Point * > & Ds2x_pts) const788 void NOMAD::TGP_Model_Search::create_Ds2x_pts
789 ( const std::vector<NOMAD::Eval_Point *> & XX , // IN
790 const NOMAD::Display & out , // IN
791 NOMAD::dd_type display_degree , // IN
792 int display_lim , // IN
793 std::vector<NOMAD::Point *> & Ds2x_pts ) const // OUT
794 {
795 if ( display_degree == NOMAD::FULL_DISPLAY )
796 out << NOMAD::open_block ( "Ds2x points construction" )
797 << std::endl;
798
799 // reset Ds2x points:
800 Ds2x_pts.clear();
801
802 // display Ds2x points directly from the NOMAD::TGP_Model object:
803 #ifdef TGP_DEBUG
804 out << NOMAD::open_block ( "expected reduction in predictive variance (Ds2x)" );
805 _model->display_Ds2x ( out , display_lim );
806 out << NOMAD::close_block() << std::endl;
807 #endif
808
809 std::set<int> pts_indexes;
810 _model->get_Ds2x_points ( pts_indexes );
811
812 if ( pts_indexes.empty() ) {
813 if ( display_degree == NOMAD::FULL_DISPLAY )
814 out << std::endl
815 << NOMAD::close_block ( "no Ds2x candidate" )
816 << std::endl;
817 return;
818 }
819
820 std::set<int>::const_iterator it , end = pts_indexes.end();
821 int i , max_index = -1 , m = _p.get_bb_nb_outputs() , n_XX = XX.size();
822
823 if ( display_degree == NOMAD::FULL_DISPLAY ) {
824 for ( it = pts_indexes.begin() ; it != end ; ++it ) {
825 if ( *it > max_index )
826 max_index = *it;
827 }
828 }
829
830 for ( it = pts_indexes.begin() , i = 0 ; it != end ; ++it , ++i ) {
831
832 Ds2x_pts.push_back ( XX[*it] );
833
834 if ( display_degree == NOMAD::FULL_DISPLAY ) {
835 out << "Ds2x candidate #";
836 out.display_int_w ( i , m );
837 out << " (XX point #";
838 out.display_int_w ( *it , max_index );
839 out << "): x=( ";
840 XX[*it]->NOMAD::Point::display ( out , " " , 6 , -1 );
841 out << " )";
842 if ( *it == n_XX - 1 )
843 out << " (rejected: candidate==incumbent)";
844 out << std::endl;
845 }
846 }
847
848 if ( display_degree == NOMAD::FULL_DISPLAY )
849 out << std::endl
850 << NOMAD::close_block ( "end of Ds2x points construction" )
851 << std::endl;
852 }
853
854 /*-------------------------------------------------------------------*/
855 /* create a list of oracle points, given by the model optimization */
856 /* (private) */
857 /*-------------------------------------------------------------------*/
create_oracle_pts(const NOMAD::Cache & cache,const NOMAD::Point & incumbent,const NOMAD::Point & delta_m,const NOMAD::Display & out,NOMAD::dd_type display_degree,int display_lim,const std::vector<NOMAD::Eval_Point * > & XX,std::vector<NOMAD::Point * > & oracle_pts,bool & stop,NOMAD::stop_type & stop_reason)858 bool NOMAD::TGP_Model_Search::create_oracle_pts
859 ( const NOMAD::Cache & cache ,
860 const NOMAD::Point & incumbent ,
861 const NOMAD::Point & delta_m ,
862 const NOMAD::Display & out ,
863 NOMAD::dd_type display_degree ,
864 int display_lim ,
865 const std::vector<NOMAD::Eval_Point *> & XX ,
866 std::vector<NOMAD::Point *> & oracle_pts ,
867 bool & stop ,
868 NOMAD::stop_type & stop_reason )
869 {
870 if ( display_degree == NOMAD::FULL_DISPLAY )
871 out << NOMAD::open_block ( "oracle points construction" )
872 << std::endl;
873
874 // reset oracle points:
875 NOMAD::TGP_Model_Search::clear_pts ( oracle_pts );
876
877 int i , n_XX = XX.size();
878
879 // starting points selection:
880 const NOMAD::Eval_Point * x0s[3];
881 x0s[0] = x0s[1] = x0s[2] = NULL;
882
883 // open display block for model predictions:
884 #ifdef TGP_DEBUG
885 out << NOMAD::open_block ( "TGP predictions (XX+ZZ)");
886 int i0 = ( display_lim > 0 ) ? n_XX - display_lim : 0;
887 if ( i0 > 0 )
888 out << "..." << std::endl;
889 #endif
890
891 NOMAD::Double f_model , h_model;
892 const NOMAD::Double & h_min = _p.get_h_min();
893 NOMAD::hnorm_type h_norm = _p.get_h_norm();
894
895 for ( i = 0 ; i < n_XX ; ++i ) {
896
897 // compute model h and f values:
898 _model->eval_hf ( XX[i]->get_bb_outputs() ,
899 h_min ,
900 h_norm ,
901 h_model ,
902 f_model );
903
904 if ( h_model.is_defined() && f_model.is_defined() ) {
905
906 XX[i]->set_f ( f_model );
907 XX[i]->set_h ( h_model );
908
909 // feasible point:
910 if ( XX[i]->is_feasible ( h_min ) ) {
911 if ( !x0s[0] || f_model < x0s[0]->get_f() )
912 x0s[0] = XX[i];
913 }
914
915 // infeasible point:
916 else {
917 if ( !x0s[1] || h_model < x0s[1]->get_h() )
918 x0s[1] = XX[i];
919
920 if ( !x0s[2] || f_model < x0s[2]->get_f() )
921 x0s[2] = XX[i];
922 }
923 }
924
925 // display model predictions:
926 #ifdef TGP_DEBUG
927 if ( i >= i0 ) {
928 out << "#";
929 out.display_int_w ( i , n_XX );
930 out << " x=(";
931 XX[i]->NOMAD::Point::display ( out , " " , 15 , -1 );
932 out << " ) m(x)=[";
933 XX[i]->get_bb_outputs().display ( out , " " , 15 , -1 );
934 out << " ]";
935
936 if ( h_model.is_defined() && f_model.is_defined() )
937 out << " hm=" << h_model << " fm=" << f_model;
938 else
939 out << " no model value";
940 out << std::endl;
941 }
942 #endif
943
944 #ifdef MODEL_STATS
945 if ( XX[i] && f_model.is_defined() && h_model.is_defined() ) {
946 XX[i]->set_mod_use ( 1 ); // 1 for model search
947 XX[i]->set_Yw ( _model->get_Yw () );
948 XX[i]->set_nY ( p );
949 XX[i]->set_mh ( h_model );
950 XX[i]->set_mf ( f_model );
951 }
952 #endif
953 }
954
955 #ifdef TGP_DEBUG
956 // close display block for model predictions:
957 {
958 std::ostringstream oss;
959 oss << "(size=" << n_XX << ")";
960 out << NOMAD::close_block ( oss.str() ) << std::endl;
961 }
962
963 // compute and display prediction errors:
964 out << NOMAD::open_block ( "prediction relative errors on X(%)" );
965 _model->display_X_errors ( out );
966 out << NOMAD::close_block() << std::endl;
967
968 #endif
969
970 if ( !x0s[0] && !x0s[1] && !x0s[2] ) {
971 if ( display_degree == NOMAD::FULL_DISPLAY )
972 out << std::endl
973 << NOMAD::close_block ( "oracle points error: no model starting point" )
974 << std::endl;
975 return false;
976 }
977
978 // display starting points:
979 if ( display_degree == NOMAD::FULL_DISPLAY ) {
980 out << std::endl
981 << NOMAD::open_block ( "model starting points" );
982
983 for ( i = 0 ; i < 3 ; ++i ) {
984
985 out << "#" << i << ": ";
986 if ( !x0s[i] )
987 out << "NULL" << std::endl;
988 else {
989
990 out << " x=(";
991 x0s[i]->NOMAD::Point::display ( out , " " , 15 , -1 );
992 out << " ) m(x)=[";
993 x0s[i]->get_bb_outputs().display ( out , " " , 15 , -1 );
994 out << " ]"
995 << " hm=" << std::setw(15) << x0s[i]->get_h()
996 << " fm=" << std::setw(15) << x0s[i]->get_f()
997 << std::endl;
998 }
999 }
1000 out << NOMAD::close_block() << std::endl;
1001 }
1002
1003 // optimize model:
1004 // ---------------
1005 NOMAD::Point * xf = NULL , * xi = NULL;
1006 NOMAD::Clock clock;
1007
1008 bool optimization_ok = optimize_model ( x0s ,
1009 out ,
1010 display_degree ,
1011 xf ,
1012 xi ,
1013 stop ,
1014 stop_reason );
1015
1016 _one_search_stats.add_optimization_time ( clock.get_CPU_time() );
1017
1018 if ( stop || !optimization_ok || ( xf == NULL && xi == NULL ) ) {
1019 std::string error_str;
1020 if ( xf == NULL && xi == NULL )
1021 error_str = "no model optimization solution";
1022 else {
1023 if ( xf ) delete xf;
1024 if ( xi ) delete xi;
1025 error_str = ( stop ) ? "algorithm stop" : "model optimization error";
1026 }
1027
1028 if ( display_degree == NOMAD::FULL_DISPLAY )
1029 out << std::endl
1030 << NOMAD::close_block ( "oracle points error: " + error_str )
1031 << std::endl;
1032 return false;
1033 }
1034
1035 // project and check xf and xi:
1036 if ( xf && !check_oracle_point ( cache ,
1037 incumbent ,
1038 delta_m ,
1039 out ,
1040 display_degree ,
1041 *xf ) ) {
1042 delete xf;
1043 xf = NULL;
1044 }
1045
1046 if ( xi && !check_oracle_point ( cache ,
1047 incumbent ,
1048 delta_m ,
1049 out ,
1050 display_degree ,
1051 *xi ) ) {
1052 delete xi;
1053 xi = NULL;
1054 }
1055
1056 // add xf and xi in the list of oracle points:
1057 if ( xf )
1058 oracle_pts.push_back ( xf );
1059
1060 if ( xi ) {
1061
1062 // check that xi != xf:
1063 if ( xf && *xf == *xi ) {
1064 delete xi;
1065 xi = NULL;
1066 }
1067
1068 if ( xi )
1069 oracle_pts.push_back ( xi );
1070 }
1071
1072 if ( display_degree == NOMAD::FULL_DISPLAY )
1073 out << std::endl
1074 << NOMAD::close_block ( "end of oracle points construction" )
1075 << std::endl;
1076
1077 return true;
1078 }
1079
1080 /*------------------------------------------------------*/
1081 /* project and accept or reject an oracle trial point */
1082 /* (private) */
1083 /*------------------------------------------------------*/
check_oracle_point(const NOMAD::Cache & cache,const NOMAD::Point & incumbent,const NOMAD::Point & delta_m,const NOMAD::Display & out,NOMAD::dd_type display_degree,NOMAD::Point & x)1084 bool NOMAD::TGP_Model_Search::check_oracle_point
1085 ( const NOMAD::Cache & cache ,
1086 const NOMAD::Point & incumbent ,
1087 const NOMAD::Point & delta_m ,
1088 const NOMAD::Display & out ,
1089 NOMAD::dd_type display_degree ,
1090 NOMAD::Point & x )
1091 {
1092 bool proj_to_mesh = _p.get_model_search_proj_to_mesh();
1093
1094 if ( display_degree == NOMAD::FULL_DISPLAY ) {
1095 out << std::endl << "oracle candidate";
1096 if ( proj_to_mesh )
1097 out << " (before projection)";
1098 out << ": ( " << x << " )" << std::endl;
1099 }
1100
1101 // projection to mesh:
1102 if ( proj_to_mesh ) {
1103 x.project_to_mesh ( incumbent , delta_m , _p.get_lb() , _p.get_ub() );
1104 if ( display_degree == NOMAD::FULL_DISPLAY )
1105 out << "oracle candidate (after projection) : ( "
1106 << x << " )" << std::endl;
1107 }
1108
1109 // compare x and incumbent coordinates:
1110 if ( x == incumbent ) {
1111 if ( display_degree == NOMAD::FULL_DISPLAY )
1112 out << "oracle candidate rejected (candidate==incumbent)" << std::endl;
1113 return false;
1114 }
1115
1116 // two evaluations points are created in order to:
1117 // 1. check if the candidate is in cache
1118 // 2. have a prediction at x and at the incumbent:
1119 int n = x.size() , m = _p.get_bb_nb_outputs();
1120
1121 NOMAD::Eval_Point * tk = new NOMAD::Eval_Point ( n , m ); // trial point
1122 tk->Point::operator = ( x );
1123
1124 // check if the point is in cache:
1125 if ( cache.find ( *tk ) ) {
1126 if ( display_degree == NOMAD::FULL_DISPLAY )
1127 out << "oracle candidate rejected (found in cache)" << std::endl;
1128 delete tk;
1129 return false;
1130 }
1131
1132 NOMAD::Eval_Point * ic = new NOMAD::Eval_Point ( n , m ); // incumbent copy
1133 ic->Point::operator = ( incumbent );
1134
1135 // model predictions (in order to accept or reject the trial point):
1136 bool pred_error = !_model->predict ( *ic , true ) || // pred_outside_bnds = true
1137 !_model->predict ( *tk , true ); // pred_outside_bnds = true
1138
1139 const NOMAD::Double & h_min = _p.get_h_min();
1140 NOMAD::hnorm_type h_norm = _p.get_h_norm();
1141
1142 NOMAD::Double h0 , f0; // model values of f and h at the center
1143 NOMAD::Double h1 , f1; // model values of f and h at the trial point
1144
1145 if ( !pred_error )
1146 _model->eval_hf ( tk->get_bb_outputs() , h_min , h_norm , h1 , f1 );
1147
1148 delete tk;
1149
1150 if ( pred_error || !h1.is_defined() || !f1.is_defined() ) {
1151
1152 if ( display_degree == NOMAD::FULL_DISPLAY ) {
1153 if ( pred_error )
1154 out << "prediction error: oracle candidate rejected";
1155 else
1156 out << "no model value (EB constraint violated?): oracle candidate rejected";
1157 out << std::endl;
1158 }
1159
1160 delete ic;
1161
1162 return false;
1163 }
1164
1165 // accept or reject the trial point:
1166 bool accept_point = false;
1167 _model->eval_hf ( ic->get_bb_outputs(), h_min, h_norm, h0, f0 );
1168 if ( !h0.is_defined() || !f0.is_defined() )
1169 accept_point = true;
1170
1171 delete ic;
1172
1173 if ( display_degree == NOMAD::FULL_DISPLAY )
1174 out << std::endl
1175 << "incumbent prediction : h0=" << h0 << " f0=" << f0 << std::endl
1176 << "trial point prediction: h1=" << h1 << " f1=" << f1 << std::endl;
1177
1178 if ( !accept_point )
1179 accept_point = (f1 < f0) || (h1 < h0);
1180
1181 if ( !accept_point ) {
1182 if ( display_degree == NOMAD::FULL_DISPLAY )
1183 out << "oracle candidate rejected" << std::endl;
1184 _one_search_stats.add_MS_rejected();
1185 return false;
1186 }
1187
1188 if ( display_degree == NOMAD::FULL_DISPLAY )
1189 out << "oracle candidate accepted" << std::endl;
1190
1191 return true;
1192 }
1193
1194 /*--------------------------------------------------------*/
1195 /* insert a trial point in the evaluator control object */
1196 /* (private) */
1197 /*--------------------------------------------------------*/
register_point(NOMAD::Point x,NOMAD::Signature & signature,const NOMAD::Point & incumbent,NOMAD::dd_type display_degree,NOMAD::Evaluator_Control & ev_control) const1198 void NOMAD::TGP_Model_Search::register_point
1199 ( NOMAD::Point x ,
1200 NOMAD::Signature & signature ,
1201 const NOMAD::Point & incumbent ,
1202 NOMAD::dd_type display_degree ,
1203 NOMAD::Evaluator_Control & ev_control ) const
1204 {
1205 int n = x.size();
1206
1207 NOMAD::Eval_Point * tk = new NOMAD::Eval_Point ( n , _p.get_bb_nb_outputs() );
1208
1209 // if the search is optimistic, a direction is computed (this
1210 // will be used in case of success in the speculative search):
1211 if ( _p.get_model_search_optimistic() ) {
1212 NOMAD::Direction dir ( n , 0.0 , NOMAD::MODEL_SEARCH_DIR );
1213 dir.Point::operator = ( x - incumbent );
1214 tk->set_direction ( &dir );
1215 }
1216
1217 tk->set_signature ( &signature );
1218 tk->Point::operator = ( x );
1219
1220 // add the new point to the list of search trial points:
1221 ev_control.add_eval_point ( tk ,
1222 display_degree ,
1223 _p.get_snap_to_bounds() ,
1224 NOMAD::Double() ,
1225 NOMAD::Double() ,
1226 NOMAD::Double() ,
1227 NOMAD::Double() );
1228 #ifdef MODEL_STATS
1229 if ( tk ) {
1230 NOMAD::Double h1 , f1;
1231 if ( _model->predict ( *tk , true ) ) // pred_outside_bnds = true
1232 _model->eval_hf ( tk->get_bb_outputs() ,
1233 _p.get_h_min() ,
1234 _p.get_h_norm() ,
1235 h1 ,
1236 f1 );
1237 if ( h1.is_defined() && f1.is_defined() ) {
1238 tk->set_mod_use ( 1 ); // 1 for model search
1239 tk->set_Yw ( _model->get_Yw() );
1240 tk->set_nY ( model.get_p () );
1241 tk->set_mh ( h1 );
1242 tk->set_mf ( f1 );
1243 }
1244 }
1245 #endif
1246 }
1247
1248 /*---------------------------------------------------------------*/
1249 /* model construction (private) */
1250 /*---------------------------------------------------------------*/
model_construction(const NOMAD::Cache & cache,const NOMAD::Point & incumbent,const NOMAD::Point & delta_m,const NOMAD::Display & out,NOMAD::dd_type display_degree,int display_lim,NOMAD::Stats & stats,bool & compute_Ds2x,std::vector<NOMAD::Eval_Point * > & XX,bool & stop,NOMAD::stop_type & stop_reason,std::string & error_str)1251 bool NOMAD::TGP_Model_Search::model_construction
1252 ( const NOMAD::Cache & cache ,
1253 const NOMAD::Point & incumbent ,
1254 const NOMAD::Point & delta_m ,
1255 const NOMAD::Display & out ,
1256 NOMAD::dd_type display_degree ,
1257 int display_lim ,
1258 NOMAD::Stats & stats ,
1259 bool & compute_Ds2x ,
1260 std::vector<NOMAD::Eval_Point *> & XX ,
1261 bool & stop ,
1262 NOMAD::stop_type & stop_reason ,
1263 std::string & error_str )
1264 {
1265 int i , n_XX , n = incumbent.size();
1266
1267 compute_Ds2x = false;
1268
1269 // reset XX:
1270 {
1271 n_XX = XX.size();
1272 for ( i = 0 ; i < n_XX ; ++i )
1273 delete XX[i];
1274 XX.clear();
1275 n_XX = 0;
1276 }
1277
1278 error_str.clear();
1279
1280 if ( stop )
1281 return false;
1282
1283 NOMAD::Clock clock;
1284
1285 // TGP model creation:
1286 if ( _model )
1287 delete _model;
1288
1289 NOMAD::TGP_mode_type tgp_mode = _p.get_model_tgp_mode();
1290
1291 _model = new NOMAD::TGP_Model ( n ,
1292 _p.get_bb_output_type() ,
1293 out ,
1294 tgp_mode );
1295
1296 // construct interpolation set (X):
1297 // --------------------------------
1298 if ( !_model->set_X ( cache ,
1299 &incumbent ,
1300 _p.get_seed() ,
1301 true ) ) { // remove_fv = true
1302
1303 if ( _model->get_nep_flag() )
1304 _one_search_stats.add_not_enough_pts();
1305
1306 stats.update_model_stats ( _one_search_stats );
1307 _all_searches_stats.update ( _one_search_stats );
1308
1309 error_str = _model->get_error_str();
1310
1311 delete _model;
1312 _model = NULL;
1313
1314 return false;
1315 }
1316
1317 // create the list of prediction points (XX)
1318 // (they are used to get starting points
1319 // and to get IMPROV and DS2X stats):
1320 // -----------------------------------
1321 set_XX ( cache ,
1322 n ,
1323 _p.get_bb_output_type().size() ,
1324 incumbent ,
1325 delta_m ,
1326 XX );
1327
1328 n_XX = XX.size();
1329
1330 // display sets X and XX:
1331 // ----------------------
1332 #ifdef TGP_DEBUG
1333
1334 // X:
1335 _model->display_X ( out , display_lim );
1336
1337 // XX (only the 10 last points are displayed):
1338 out << NOMAD::open_block ( "prediction points (XX)");
1339 int i0 = ( display_lim > 0 ) ? n_XX - display_lim : 0;
1340 if ( i0 > 0 )
1341 out << "..." << std::endl;
1342 else if ( i0 < 0 )
1343 i0 = 0;
1344 for ( i = i0 ; i < n_XX ; ++i ) {
1345 out << "#";
1346 out.display_int_w ( i , n_XX );
1347 out << " x=(";
1348 XX[i]->NOMAD::Point::display ( out , " " , 15 , -1 );
1349 out << " )" << std::endl;
1350 }
1351 std::ostringstream oss;
1352 oss << "(size=" << n_XX << ")";
1353 out << NOMAD::close_block ( oss.str() ) << std::endl;
1354 #endif
1355
1356 // model construction:
1357 // -------------------
1358 int p = _model->get_p(); // number of points in X
1359
1360 if ( display_degree == NOMAD::FULL_DISPLAY ) {
1361 out << "TGP model construction (p="
1362 << p << ") ...";
1363 out.flush();
1364 }
1365
1366 // decide compute_Ds2x flag:
1367 compute_Ds2x = true;
1368 if ( tgp_mode == NOMAD::TGP_PRECISE && n_XX > 100 )
1369 compute_Ds2x = false;
1370
1371 if ( !_model->compute ( XX ,
1372 compute_Ds2x ,
1373 true , // compute_improv = true
1374 false ) ) { // pred_outside_bnds = false
1375
1376 _one_search_stats.add_construction_error();
1377
1378 if ( display_degree == NOMAD::FULL_DISPLAY )
1379 out << "... error" << std::endl;
1380
1381 // delete XX points:
1382 for ( i = 0 ; i < n_XX ; ++i )
1383 delete XX[i];
1384 XX.clear();
1385 n_XX = 0;
1386
1387 stats.update_model_stats ( _one_search_stats );
1388 _all_searches_stats.update ( _one_search_stats );
1389
1390 error_str = _model->get_error_str();
1391
1392 delete _model;
1393 _model = NULL;
1394
1395 // check if ctrl-c has been pressed:
1396 if ( NOMAD::TGP_Output_Model::get_force_quit() ) {
1397 stop = true;
1398 stop_reason = NOMAD::CTRL_C;
1399 }
1400
1401 return false;
1402 }
1403
1404 if ( display_degree == NOMAD::FULL_DISPLAY )
1405 out << "... OK" << std::endl << std::endl;
1406
1407 // update model stats:
1408 _one_search_stats.add_construction_time ( clock.get_CPU_time() );
1409 _one_search_stats.update_nY ( p );
1410 _one_search_stats.add_nb_truth();
1411 _one_search_stats.add_nb_TGP();
1412
1413 return true;
1414 }
1415
1416 /*---------------------------------------------------------------*/
1417 /* optimize a model (private) */
1418 /*---------------------------------------------------------------*/
optimize_model(const NOMAD::Eval_Point * x0s[3],const NOMAD::Display & out,NOMAD::dd_type display_degree,NOMAD::Point * & xf,NOMAD::Point * & xi,bool & stop,NOMAD::stop_type & stop_reason)1419 bool NOMAD::TGP_Model_Search::optimize_model
1420 ( const NOMAD::Eval_Point * x0s[3] ,
1421 const NOMAD::Display & out ,
1422 NOMAD::dd_type display_degree ,
1423 NOMAD::Point *& xf ,
1424 NOMAD::Point *& xi ,
1425 bool & stop ,
1426 NOMAD::stop_type & stop_reason )
1427 {
1428 // reset xf and xi:
1429 if ( xf ) delete xf; xf = NULL;
1430 if ( xi ) delete xi; xi = NULL;
1431
1432 if ( !_model ) {
1433 if ( display_degree == NOMAD::FULL_DISPLAY )
1434 out << std::endl << "model optimization error (no model)"
1435 << std::endl;
1436 return false;
1437 }
1438
1439 std::string error_str;
1440 bool error = false;
1441 int i , n = _model->get_n0();
1442
1443 // model bounds:
1444 const NOMAD::Point & lb = _model->get_lb();
1445 const NOMAD::Point & ub = _model->get_ub();
1446
1447 // initial displays:
1448 if ( display_degree == NOMAD::FULL_DISPLAY )
1449 out << std::endl << NOMAD::open_block ( "model optimization" );
1450
1451 // parameters creation:
1452 NOMAD::Parameters model_param ( out );
1453
1454
1455 // number of variables:
1456 model_param.set_DIMENSION ( n );
1457
1458 // blackbox outputs:
1459 model_param.set_BB_OUTPUT_TYPE ( _p.get_bb_output_type() );
1460
1461 // barrier parameters:
1462 model_param.set_H_MIN ( _p.get_h_min () );
1463 model_param.set_H_NORM ( _p.get_h_norm() );
1464
1465 // starting points:
1466 for ( i = 0 ; i < 3 ; ++i )
1467 if ( x0s[i] )
1468 model_param.set_X0 ( *x0s[i] );
1469
1470 // fixed variables:
1471 for ( i = 0 ; i < n ; ++i )
1472 if ( lb[i] == ub[i] || _p.variable_is_fixed(i) )
1473 model_param.set_FIXED_VARIABLE(i);
1474
1475 // no model search and no model ordering:
1476 model_param.set_MODEL_SEARCH ( false );
1477 model_param.set_MODEL_EVAL_SORT ( false );
1478
1479 // display:
1480 model_param.set_DISPLAY_DEGREE ( NOMAD::NO_DISPLAY );
1481 if ( display_degree == NOMAD::FULL_DISPLAY ) {
1482
1483 model_param.set_DISPLAY_DEGREE ( NOMAD::NORMAL_DISPLAY );
1484 // model_param.set_DISPLAY_DEGREE ( NOMAD::FULL_DISPLAY );
1485
1486 if ( n <= 5 )
1487 model_param.set_DISPLAY_STATS ( "bbe mesh_index ( %14.12gsol ) %14.12gobj" );
1488 else if ( n <= 10 )
1489 model_param.set_DISPLAY_STATS ( "bbe mesh_index ( sol ) obj" );
1490 else
1491 model_param.set_DISPLAY_STATS ( "bbe obj" );
1492 }
1493 // mesh: use isotropic mesh
1494 model_param.set_ANISOTROPIC_MESH ( false );
1495 model_param.set_MESH_UPDATE_BASIS ( 4.0 );
1496 model_param.set_MESH_COARSENING_EXPONENT ( 1 );
1497 model_param.set_MESH_REFINING_EXPONENT ( -1 );
1498 model_param.set_INITIAL_MESH_INDEX ( 0 );
1499
1500
1501 // maximum number of evaluations (2000 or 10000):
1502 model_param.set_MAX_BB_EVAL
1503 ( ( _p.get_model_tgp_mode() == NOMAD::TGP_PRECISE ) ? 10000 : 2000 );
1504
1505 model_param.set_SNAP_TO_BOUNDS ( true );
1506
1507 // disable user calls:
1508 model_param.set_USER_CALLS_ENABLED ( false );
1509
1510 // set flags:
1511 bool flag_check_bimads , flag_reset_mesh , flag_reset_barriers , flag_p1_active;
1512 NOMAD::Mads::get_flags ( flag_check_bimads ,
1513 flag_reset_mesh ,
1514 flag_reset_barriers ,
1515 flag_p1_active );
1516
1517 NOMAD::Mads::set_flag_check_bimads ( true );
1518 NOMAD::Mads::set_flag_reset_mesh ( true );
1519 NOMAD::Mads::set_flag_reset_barriers ( true );
1520 NOMAD::Mads::set_flag_p1_active ( false );
1521
1522 // bounds:
1523 model_param.set_LOWER_BOUND ( lb );
1524 model_param.set_UPPER_BOUND ( ub );
1525
1526 try {
1527
1528 // parameters validation:
1529 model_param.check();
1530
1531 // out << "TGP PARAMETERS:" << std::endl << model_param << std::endl;
1532
1533 // model evaluator creation:
1534 NOMAD::TGP_Model_Evaluator ev ( model_param , *_model );
1535
1536 // algorithm creation and execution:
1537 NOMAD::Mads mads ( model_param , &ev );
1538 NOMAD::stop_type st = mads.run();
1539
1540 // check the stopping criterion:
1541 if ( st == NOMAD::CTRL_C || st == NOMAD::MAX_CACHE_MEMORY_REACHED ) {
1542 std::ostringstream oss;
1543 oss << "model optimization: " << st;
1544 error_str = oss.str();
1545 error = true;
1546 stop = true;
1547 stop_reason = st;
1548 }
1549
1550 else if ( st == NOMAD::MAX_BB_EVAL_REACHED )
1551 _one_search_stats.add_MS_max_bbe();
1552
1553 // update the stats on the number of model evaluations:
1554 _one_search_stats.update_MS_model_opt ( mads.get_stats().get_bb_eval() );
1555
1556 // get the solution(s):
1557 const NOMAD::Eval_Point * best_feas = mads.get_best_feasible ();
1558 const NOMAD::Eval_Point * best_infeas = mads.get_best_infeasible();
1559
1560 if ( best_feas )
1561 xf = new NOMAD::Point ( *best_feas );
1562 if ( best_infeas )
1563 xi = new NOMAD::Point ( *best_infeas );
1564
1565 if ( !xf && !xi ) {
1566 error = true;
1567 error_str = "optimization error: no solution";
1568 }
1569 }
1570 catch ( std::exception & e ) {
1571 error = true;
1572 error_str = std::string ( "optimization error: " ) + e.what();
1573 }
1574
1575 // reset flags:
1576 NOMAD::Mads::set_flag_check_bimads ( flag_check_bimads );
1577 NOMAD::Mads::set_flag_reset_mesh ( flag_reset_mesh );
1578 NOMAD::Mads::set_flag_reset_barriers ( flag_reset_barriers );
1579 NOMAD::Mads::set_flag_p1_active ( flag_p1_active );
1580
1581 // close display block:
1582 if ( display_degree == NOMAD::FULL_DISPLAY ) {
1583 if ( error )
1584 out.close_block ( error_str );
1585 else
1586 out.close_block();
1587 }
1588
1589 return !error;
1590 }
1591
1592 /*---------------------------------------------------------*/
1593 /* display the prediction error for the evaluated points */
1594 /* (private) */
1595 /*---------------------------------------------------------*/
predict(const NOMAD::Point & x,NOMAD::Double & h,NOMAD::Double & f) const1596 bool NOMAD::TGP_Model_Search::predict ( const NOMAD::Point & x ,
1597 NOMAD::Double & h ,
1598 NOMAD::Double & f ) const
1599 {
1600 h.clear();
1601 f.clear();
1602
1603 if ( !_model )
1604 return false;
1605
1606 NOMAD::Eval_Point y ( x.size() , _p.get_bb_nb_outputs() );
1607 y.NOMAD::Point::operator = ( x );
1608
1609 if ( _model->predict ( y , true ) ) // pred_outside_bnds = true
1610 _model->eval_hf ( y.get_bb_outputs() , _p.get_h_min() , _p.get_h_norm() , h , f );
1611
1612 return ( h.is_defined() && f.is_defined() );
1613 }
1614
1615 /*---------------------------------------------------------*/
1616 /* display the prediction error for the evaluated points */
1617 /* (private) */
1618 /*---------------------------------------------------------*/
display_eval_pred_errors(const std::list<const NOMAD::Eval_Point * > & evaluated_pts,const NOMAD::Display & out)1619 void NOMAD::TGP_Model_Search::display_eval_pred_errors
1620 ( const std::list<const NOMAD::Eval_Point *> & evaluated_pts ,
1621 const NOMAD::Display & out )
1622 {
1623 if ( !_model )
1624 return;
1625
1626 int i , j = 0;
1627 int nb_pts = evaluated_pts.size() ,
1628 n = _model->get_n0() ,
1629 m = _p.get_bb_nb_outputs();
1630 const NOMAD::Double & h_min = _p.get_h_min();
1631 NOMAD::hnorm_type h_norm = _p.get_h_norm();
1632 NOMAD::Double h , f , hm , fm , err;
1633 NOMAD::Eval_Point x ( n , m );
1634
1635 out << std::endl << NOMAD::open_block ( "evaluated points" );
1636 std::list<const NOMAD::Eval_Point *>::const_iterator it , end = evaluated_pts.end();
1637 for ( it = evaluated_pts.begin() ; it != end ; ++it ) {
1638
1639 if ( !(*it) )
1640 continue;
1641
1642 h = (*it)->get_h();
1643 f = (*it)->get_f();
1644
1645 for ( i = 0 ; i < n ; ++i )
1646 x[i] = (**it)[i];
1647
1648 _model->predict ( x , true ); // pred_outside_bnds = true
1649 _model->eval_hf ( x.get_bb_outputs(), h_min , h_norm , hm , fm );
1650
1651 out << "#";
1652 out.display_int_w ( j++ , nb_pts );
1653
1654 out << " x=(";
1655 if ( n < 5 )
1656 (*it)->NOMAD::Point::display ( out , " " , 6 , -1 );
1657 else
1658 (*it)->NOMAD::Point::display ( out , " " , -1 , 20 );
1659 out << ")";
1660 if ( _p.has_constraints() ) {
1661 err = (h.is_defined() && hm.is_defined()) ? h.rel_err(hm) * 100.0 : 100.0;
1662 out << " [h=";
1663 h.display ( out , "%11.3g" );
1664 out << " hm=";
1665 hm.display ( out , "%11.3g" );
1666 out << " err_h=";
1667 err.display ( out , "%.2f" );
1668 out << "%]";
1669 }
1670 err = (f.is_defined() && fm.is_defined()) ? f.rel_err(fm) * 100.0 : 100.0;
1671 out << " [f=";
1672 f.display ( out , "%11.3g" );
1673 out << " fm=";
1674 fm.display ( out , "%11.3g" );
1675 out << " err_f=";
1676 err.display ( out , "%.2f" );
1677 out << "%]" << std::endl;
1678 }
1679
1680 out.close_block();
1681 }
1682
1683 #endif
1684