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.cpp
38 \brief TGP (Bayesian treed Gaussian process) model for all outputs (implementation)
39 \author Sebastien Le Digabel
40 \date 2011-02-07
41 \see TGP_Model.hpp
42 */
43
44 #ifndef USE_TGP
45
46 int TGP_DUMMY; // avoids that TGP_Model.o has no symbols with ranlib
47
48 #else
49
50 #include "TGP_Model.hpp"
51
52 /*-----------------------------------*/
53 /* constructor 1/2 */
54 /*-----------------------------------*/
TGP_Model(int n0,const std::vector<NOMAD::bb_output_type> & bbot,const NOMAD::Display & out,NOMAD::TGP_mode_type mode)55 NOMAD::TGP_Model::TGP_Model ( int n0 ,
56 const std::vector<NOMAD::bb_output_type> & bbot ,
57 const NOMAD::Display & out ,
58 NOMAD::TGP_mode_type mode )
59 : _out ( out ) ,
60 _bbot ( bbot ) ,
61 _p ( 0 ) ,
62 _n0 ( n0 ) ,
63 _n ( 0 ) ,
64 _n_XX ( 0 ) ,
65 _av_index ( NULL ) ,
66 _fv_index ( NULL ) ,
67 _X ( NULL ) ,
68 _XX ( NULL ) ,
69 _Ds2x ( NULL ) ,
70 _improv ( NULL ) ,
71 _model_computed ( false ) ,
72 _nep_flag ( false ) ,
73 _tgp_mode ( mode ) ,
74 _tgp_models ( NULL ) ,
75 _tgp_rect ( NULL ) ,
76 _tgp_linburn ( true ) ,
77 _usr_pmax ( -1 ) ,
78 _tgp_verb ( false )
79 {
80 _usr_BTE[0] = _usr_BTE[1] = _usr_BTE[2] = -1;
81
82 // user mode: this is not the good constructor to call:
83 if ( mode == NOMAD::TGP_USER )
84 throw NOMAD::Exception ( "TGP_Model.cpp" , __LINE__ ,
85 "this constructor only accepts fast or precise TGP modes" );
86
87 _error_str = "NOMAD::TGP_Model::set_X() has not been called";
88 }
89
90 /*--------------------------------------------------------*/
91 /* constructor 2/2 (with user BTE and pmax and no bbot) */
92 /*--------------------------------------------------------*/
TGP_Model(int n0,const NOMAD::Display & out,int BTE[3],bool linburn,int pmax,bool verb)93 NOMAD::TGP_Model::TGP_Model ( int n0 ,
94 const NOMAD::Display & out ,
95 int BTE[3] ,
96 bool linburn ,
97 int pmax ,
98 bool verb )
99 : _out ( out ) ,
100 _p ( 0 ) ,
101 _n0 ( n0 ) ,
102 _n ( 0 ) ,
103 _n_XX ( 0 ) ,
104 _av_index ( NULL ) ,
105 _fv_index ( NULL ) ,
106 _X ( NULL ) ,
107 _XX ( NULL ) ,
108 _Ds2x ( NULL ) ,
109 _improv ( NULL ) ,
110 _model_computed ( false ) ,
111 _nep_flag ( false ) ,
112 _tgp_mode ( NOMAD::TGP_USER ) ,
113 _tgp_models ( NULL ) ,
114 _tgp_rect ( NULL ) ,
115 _tgp_linburn ( linburn ) ,
116 _usr_pmax ( pmax ) ,
117 _tgp_verb ( verb )
118 {
119
120 // user BTE parameters:
121 _usr_BTE[0] = BTE[0];
122 _usr_BTE[1] = BTE[1];
123 _usr_BTE[2] = BTE[2];
124
125 // default bbot: just one objective output:
126 _bbot.push_back ( NOMAD::OBJ );
127
128 _error_str = "NOMAD::TGP_Model::set_X() has not been called";
129 }
130
131 /*--------------------------------------------*/
132 /* destructor */
133 /*--------------------------------------------*/
~TGP_Model(void)134 NOMAD::TGP_Model::~TGP_Model ( void )
135 {
136 clear();
137
138 if ( _Ds2x ) {
139 int i , m = _bbot.size();
140 for ( i = 0 ; i < m ; ++i )
141 if ( _Ds2x[i] )
142 delete [] _Ds2x[i];
143 delete [] _Ds2x;
144 }
145
146 if ( _improv )
147 delete [] _improv;
148 }
149
150 /*--------------------------------------------*/
151 /* clear memory (private) */
152 /*--------------------------------------------*/
153 /* members that depend only on constructor */
154 /* arguments are not reseted here */
155 /*--------------------------------------------*/
clear(void)156 void NOMAD::TGP_Model::clear ( void )
157 {
158 int i;
159
160 if ( _av_index ) {
161 delete [] _av_index;
162 _av_index = NULL;
163 }
164
165 if ( _fv_index ) {
166 delete [] _fv_index;
167 _fv_index = NULL;
168 }
169
170 if ( _X ) {
171 for ( i = 0 ; i < _p ; ++i )
172 delete [] _X[i];
173 delete [] _X;
174 _X = NULL;
175 }
176
177 if ( _XX ) {
178 for ( i = 0 ; i < _n_XX ; ++i )
179 delete [] _XX[i];
180 delete [] _XX;
181 _XX = NULL;
182 }
183
184 if ( _tgp_models ) {
185 int m = _bbot.size();
186 for ( i = 0 ; i < m ; ++i )
187 if ( _tgp_models[i] ) {
188 delete _tgp_models[i];
189 }
190 delete [] _tgp_models;
191 _tgp_models = NULL;
192 }
193
194 if ( _tgp_rect ) {
195 delete_matrix ( _tgp_rect );
196 _tgp_rect = NULL;
197 }
198
199 _error_str.clear();
200 _lb.clear();
201 _ub.clear();
202 _fv.clear();
203
204 _p = _n = _n_XX = 0;
205
206 _model_computed = false;
207 }
208
209 /*--------------------------------------------*/
210 /* set the interpolation set X (1/2) */
211 /*--------------------------------------------*/
set_X(const NOMAD::Cache & cache,const NOMAD::Point * center,int seed,bool remove_fv)212 bool NOMAD::TGP_Model::set_X ( const NOMAD::Cache & cache , // only truth evals
213 const NOMAD::Point * center , // the incumbent
214 int seed ,
215 bool remove_fv )
216 {
217 std::vector<const NOMAD::Eval_Point *> X;
218
219 const NOMAD::Eval_Point * cur = cache.begin();
220 while ( cur ) {
221 X.push_back ( cur );
222 cur = cache.next();
223 }
224
225 return set_X ( X , center , seed , remove_fv );
226 }
227
228 /*--------------------------------------------*/
229 /* set the interpolation set X (2/2) */
230 /*--------------------------------------------*/
set_X(const std::vector<const NOMAD::Eval_Point * > & X,const NOMAD::Point * center,int seed,bool remove_fv)231 bool NOMAD::TGP_Model::set_X
232 ( const std::vector<const NOMAD::Eval_Point *> & X ,
233 const NOMAD::Point * center ,
234 int seed ,
235 bool remove_fv )
236 {
237 clear();
238 _nep_flag = false;
239
240 if ( X.size() <= 1 ) {
241 _error_str = ( X.empty() ) ?
242 "interpolation set X is empty" : "only one point in interpolation set X";
243 _nep_flag = true;
244 return false;
245 }
246
247 // filter and sort X:
248 std::list<const NOMAD::Eval_Point *> filtered_X;
249 _p = filter_and_sort_X ( X , center , filtered_X );
250
251 if ( _p == 0 ) {
252 _error_str = "no valid interpolation point";
253 _nep_flag = true;
254 return false;
255 }
256
257 int pmin , pmax;
258 bool valid = false;
259
260 while ( !valid ) {
261
262 _n = check_fixed_variables ( filtered_X , remove_fv );
263
264 if ( _n <= 1 ) {
265 _error_str = ( _n == 0 ) ?
266 "only fixed variables" : "interpolation set of dimension 1";
267 return false;
268 }
269
270 // get the limits pmin <= p <= pmax:
271 if ( !get_p_limits ( _n , pmin , pmax ) )
272 return false;
273
274 if ( _p < pmin ) {
275 std::ostringstream oss;
276 oss << "not enough interpolation points ("
277 << _p << "<=" << pmin << ")";
278 _error_str = oss.str();
279 _nep_flag = true;
280 return false;
281 }
282
283 // reduce the number of interpolation points:
284 if ( _p > pmax ) {
285 filtered_X.resize ( pmax );
286 _p = pmax;
287 valid = false;
288 }
289
290 // the interpolation set has a valid size:
291 else
292 valid = true;
293 }
294
295 // display limits on p:
296 #ifdef TGP_DEBUG
297 _out << std::endl
298 << "number of interpolation points: "
299 << pmin << " <= p=" << _p << " <= " << pmax
300 << std::endl;
301 #endif
302
303 // create interpolation matrix X:
304 int i = 0 , j;
305 _X = new double * [_p];
306
307 std::list<const NOMAD::Eval_Point *>::const_iterator
308 it , end = filtered_X.end();
309 for ( it = filtered_X.begin() ; it != end ; ++it , ++i ) {
310 _X[i] = new double[_n];
311 for ( j = 0 ; j < _n ; ++j )
312 _X[i][j] = (**it)[_fv_index[j]].value();
313 }
314
315 // check that X is of full rank:
316 if ( !NOMAD::TGP_Model::check_full_rank ( _X , _p , _n ) ) {
317 _error_str = "X is not of full rank";
318 return false;
319 }
320
321 // construct individual TGP models (one for each output):
322 int m = _bbot.size();
323 _tgp_models = new NOMAD::TGP_Output_Model * [m];
324 for ( i = 0 ; i < m ; ++i ) {
325 _tgp_models[i] = NULL;
326 if ( _bbot[i] == NOMAD::OBJ || NOMAD::bbot_is_constraint(_bbot[i]) )
327 _tgp_models[i] = new NOMAD::TGP_Output_Model ( filtered_X ,
328 i ,
329 seed ,
330 _out );
331 }
332
333 return _error_str.empty();
334 }
335
336 /*----------------------------------------------------*/
337 /* get the limits pmin <= p <= pmax (private) */
338 /*----------------------------------------------------*/
get_p_limits(int n,int & pmin,int & pmax)339 bool NOMAD::TGP_Model::get_p_limits ( int n , int & pmin , int & pmax )
340 {
341 if ( n <= 0 ) {
342 _error_str = "bad value of n for the computation of pmin and pmax";
343 return false;
344 }
345
346 // TGP fast mode:
347 // --------------
348 //
349 // example of some values for pmin and pmax:
350 //
351 // n pmin pmax
352 //
353 // 2 7 30
354 // 3 8 30
355 // 4 9 30
356 // 5 10 30
357 // 6 11 32
358 // 7 12 35
359 // 8 13 37
360 // 9 14 40
361 // 10 15 42
362 // 20 25 60
363 // 50 55 94
364 // 55 60 99
365 // 56 61 100
366 // 57 62 100
367 // 95 100 100
368 // 96 101 101
369 // 100 105 105
370
371 if ( _tgp_mode == NOMAD::TGP_FAST ) {
372 pmin = n+5;
373
374 pmax = static_cast<int> ( floor ( sqrt(180.0*n) ) );
375
376 if ( pmax > 100 )
377 pmax = 100;
378
379 if ( pmax < 30 )
380 pmax = 30;
381
382 if ( pmax < pmin )
383 pmax = pmin;
384
385 return true;
386 }
387
388 // TGP precise mode: pmin=n+1 and pmax=100
389 // -----------------
390 if ( _tgp_mode == NOMAD::TGP_PRECISE ) {
391 pmin = n+1;
392 pmax = 100;
393 return true;
394 }
395
396 // TGP user mode: pmin=n+1 and pmax is user decided with constructor 2/2
397 // --------------
398 if ( _tgp_mode == NOMAD::TGP_USER ) {
399 pmin = n+1;
400 pmax = _usr_pmax;
401 if ( pmax < pmin ) {
402 std::ostringstream oss;
403 oss << "user pmax value must be > " << n;
404 _error_str = oss.str();
405 return false;
406 }
407 return true;
408 }
409
410 return true;
411 }
412
413 /*----------------------------------------------*/
414 /* get the TGP BTE parameters (private) */
415 /*----------------------------------------------*/
get_BTE(int BTE[3])416 bool NOMAD::TGP_Model::get_BTE ( int BTE[3] )
417 {
418 // fast TGP mode:
419 // --------------
420 // . if p=pmin, then B=2000
421 // . if p=pmax, then B=1000
422 // . if pmin<p<pmax, a linear relation is used
423 // . B is rounded so that 10 divides it
424 // . T=3B
425 // . E=10 (E divides T-B)
426 if ( _tgp_mode == NOMAD::TGP_FAST ) {
427
428 int pmin , pmax;
429 get_p_limits ( _n , pmin , pmax );
430
431 double a = 1000.0 / (pmin-pmax);
432 double b = 2000 - a * pmin;
433
434 BTE[0] = static_cast<int> ( 10.0 * floor ( (a * _p + b) / 10.0 ) );
435 BTE[1] = 3*BTE[0];
436 BTE[2] = 10;
437 }
438
439 // precise TGP mode:
440 // -----------------
441 else if ( _tgp_mode == NOMAD::TGP_PRECISE ) {
442 BTE[0] = 2000;
443 BTE[1] = 7000;
444 BTE[2] = 2;
445 }
446
447 // user mode:
448 // ----------
449 else {
450 BTE[0] = _usr_BTE[0];
451 BTE[1] = _usr_BTE[1];
452 BTE[2] = _usr_BTE[2];
453
454 // check the user BTE parameters:
455 if ( BTE[0] <= 0 ||
456 BTE[1] <= 0 ||
457 BTE[2] <= 0 ||
458 BTE[1] <= BTE[0] ) {
459 _error_str = "error with user BTE";
460 return false;
461 }
462 }
463
464 return true;
465 }
466
467 /*--------------------------------------------------------------*/
468 /* compute the ranges, check the fixed variables, set the */
469 /* different indexes, and return the number of free variables */
470 /* (private) */
471 /*--------------------------------------------------------------*/
check_fixed_variables(const std::list<const NOMAD::Eval_Point * > & X,bool remove_fv)472 int NOMAD::TGP_Model::check_fixed_variables
473 ( const std::list<const NOMAD::Eval_Point *> & X ,
474 bool remove_fv )
475 {
476
477 int i;
478
479 if ( !_av_index ) {
480 _av_index = new int [_n0];
481 for ( i = 0 ; i < _n0 ; ++i )
482 _av_index[i] = -1;
483 }
484
485 _lb = NOMAD::Point ( _n0 );
486 _ub = NOMAD::Point ( _n0 );
487 _fv = NOMAD::Point ( _n0 );
488
489 // compute ranges:
490 std::list<const NOMAD::Eval_Point *>::const_iterator it , end = X.end();
491 for ( it = X.begin() ; it != end ; ++it ) {
492 for ( i = 0 ; i < _n0 ; ++i ) {
493 if ( !_lb[i].is_defined() || (**it)[i] < _lb[i] )
494 _lb[i] = (**it)[i];
495 if ( !_ub[i].is_defined() || (**it)[i] > _ub[i] )
496 _ub[i] = (**it)[i];
497 }
498 }
499
500 // compute n (number of free variables), the fixed variables, and the indexes:
501 int n = 0;
502
503 std::vector<int> fv_index_tmp;
504 for ( i = 0 ; i < _n0 ; ++i ) {
505 if ( remove_fv && _lb[i] == _ub[i] ) {
506 _fv [i] = _lb[i];
507 _av_index[i] = -1;
508 }
509 else {
510 fv_index_tmp.push_back ( i );
511 _av_index[i] = n++;
512 }
513 }
514
515 // complete the fixed var index construction:
516 if ( _fv_index )
517 delete [] _fv_index;
518 _fv_index = NULL;
519
520 if ( n > 0 ) {
521 _fv_index = new int[n];
522 for ( i = 0 ; i < n ; ++i )
523 _fv_index[i] = fv_index_tmp[i];
524 }
525
526 return n;
527 }
528
529 /*--------------------------------------------------------*/
530 /* filter and sort an interpolation set X (private) */
531 /* . the points are sorted relatively to the distance */
532 /* from the center */
533 /* . if the center is not defined (equal to NULL), */
534 /* an alternate center is constructed */
535 /*--------------------------------------------------------*/
filter_and_sort_X(const std::vector<const NOMAD::Eval_Point * > & X,const NOMAD::Point * center,std::list<const NOMAD::Eval_Point * > & filtered_X) const536 int NOMAD::TGP_Model::filter_and_sort_X
537 ( const std::vector<const NOMAD::Eval_Point *> & X ,
538 const NOMAD::Point * center ,
539 std::list<const NOMAD::Eval_Point *> & filtered_X ) const
540 {
541 NOMAD::Point alt_center;
542 const NOMAD::Eval_Point * cur = NULL;
543 int p0 = X.size() , i;
544
545 // alternate center if center==NULL:
546 if ( !center ) {
547 int j;
548 NOMAD::Point lb(_n0) , ub(_n0);
549 for ( i = 0 ; i < p0 ; ++i ) {
550 cur = X[i];
551 if ( test_interpolation_point ( cur ) ) {
552 for ( j = 0 ; j < _n0 ; ++j ) {
553 if ( !lb[j].is_defined() || (*cur)[j] < lb[j] )
554 lb[j] = (*cur)[j];
555 if ( !ub[j].is_defined() || (*cur)[j] > ub[j] )
556 ub[j] = (*cur)[j];
557 }
558 }
559 }
560 alt_center = NOMAD::Point(_n0);
561 for ( j = 0 ; j < _n0 ; ++j )
562 alt_center[j] = ( lb[j] + ub[j] ) / 2.0;
563 }
564
565 // X_tmp is used to sort the points:
566 std::multiset<NOMAD::Model_Sorted_Point> tmp_X;
567
568 for ( i = 0 ; i < p0 ; ++i ) {
569
570 cur = X[i];
571
572 // test if the interpolation point is valid for interpolation:
573 if ( test_interpolation_point ( cur ) ) {
574
575 NOMAD::Model_Sorted_Point sorted_pt
576 ( &NOMAD::Cache::get_modifiable_point (*cur) ,
577 (center) ? *center : alt_center );
578
579 tmp_X.insert ( sorted_pt );
580 }
581 }
582
583 // copy the set X_tmp to filtered_X:
584 std::multiset<NOMAD::Model_Sorted_Point>::const_iterator it , end = tmp_X.end();
585 for ( it = tmp_X.begin() ; it != end ; ++it )
586 filtered_X.push_back ( static_cast<NOMAD::Eval_Point *> ( it->get_point() ) );
587
588 return filtered_X.size();
589 }
590
591 /*-----------------------------------------------------*/
592 /* tests to check if an interpolation point is valid */
593 /* (private) */
594 /*-----------------------------------------------------*/
test_interpolation_point(const NOMAD::Eval_Point * x) const595 bool NOMAD::TGP_Model::test_interpolation_point ( const NOMAD::Eval_Point * x ) const
596 {
597 if ( !x || x->size() != _n0 || !x->is_eval_ok() )
598 return false;
599
600 int m = _bbot.size();
601 const NOMAD::Point & bbo = x->get_bb_outputs();
602
603 if ( bbo.size() != m )
604 return false;
605
606 for ( int j = 0 ; j < m ; ++j )
607 if ( ( _bbot[j] == NOMAD::OBJ || NOMAD::bbot_is_constraint(_bbot[j]) ) &&
608 ( !bbo[j].is_defined() || bbo[j].abs() > NOMAD::MODEL_MAX_OUTPUT ) )
609 return false;
610
611 // point is valid:
612 return true;
613 }
614
615 /*--------------------------------------------*/
616 /* compute the models (one for each output) */
617 /*--------------------------------------------*/
compute(std::vector<NOMAD::Eval_Point * > & XX_pts,bool compute_Ds2x,bool compute_improv,bool pred_outside_bnds)618 bool NOMAD::TGP_Model::compute
619 ( std::vector<NOMAD::Eval_Point *> & XX_pts , // IN/OUT
620 bool compute_Ds2x , // IN
621 bool compute_improv , // IN
622 bool pred_outside_bnds ) // IN
623 {
624 _model_computed = false;
625
626 if ( !_error_str.empty() )
627 return false;
628
629 int i , j , index_obj = -1 , n_XX0 = XX_pts.size() , m = _bbot.size();
630
631 // check bbot: there must be exactly one objective:
632 for ( i = 0 ; i < m ; ++i ) {
633 if ( _bbot[i] == NOMAD::OBJ ) {
634 if ( index_obj < 0 )
635 index_obj = i;
636 else {
637 _error_str = "more than one objective";
638 return false;
639 }
640 }
641 }
642 if ( index_obj < 0 ) {
643 _error_str = "no objective";
644 return false;
645 }
646
647 // check n_XX0:
648 if ( n_XX0 == 0 ) {
649 _error_str = "no user-provided prediction point";
650 return false;
651 }
652
653 // reset XX_pts outputs:
654 for ( i = 0 ; i < n_XX0 ; ++i ) {
655 for ( j = 0 ; j < m ; ++j )
656 XX_pts[i]->set_bb_output ( j , NOMAD::Double() );
657 XX_pts[i]->set_eval_status ( NOMAD::EVAL_FAIL );
658 }
659
660 // create the XX matrix (prediction points):
661 std::vector<NOMAD::Eval_Point *> XX_filtered_pts;
662 NOMAD::Eval_Point * cur2;
663
664 // the list of points has to be filtered:
665 NOMAD::Double tmp;
666 bool chk;
667
668 for ( i = 0 ; i < n_XX0 ; ++i ) {
669 if ( XX_pts[i]->size() == _n0 ) {
670 cur2 = XX_pts[i];
671 chk = true;
672 for ( j = 0 ; j < _n0 ; ++j ) {
673 tmp = (*cur2)[j];
674 if ( !pred_outside_bnds && ( tmp < _lb[j] || tmp > _ub[j] ) ) {
675 chk = false;
676 break;
677 }
678 }
679 if ( chk )
680 XX_filtered_pts.push_back ( cur2 );
681 }
682 }
683
684 if ( _XX ) {
685 for ( i = 0 ; i < _n_XX ; ++i )
686 delete [] _XX[i];
687 delete [] _XX;
688 }
689
690 _n_XX = XX_filtered_pts.size();
691
692 if ( _n_XX == 0 ) {
693 _error_str = "no prediction point after filtering";
694 return false;
695 }
696
697 _XX = new double * [_n_XX];
698 for ( i = 0 ; i < _n_XX ; ++i ) {
699 _XX[i] = new double[_n];
700 for ( j = 0 ; j < _n ; ++j )
701 _XX[i][j] = (*XX_filtered_pts[i])[_fv_index[j]].value();
702 }
703
704 // Xsplit: X+XX: size = nsplit x n:
705 int nsplit = _p + _n_XX;
706 double ** Xsplit = new double * [nsplit];
707
708 for ( i = 0 ; i < _p ; ++i ) {
709 Xsplit[i] = new double [_n];
710 for ( j = 0 ; j < _n ; ++j )
711 Xsplit[i][j] = _X[i][j];
712 }
713
714 for ( i = _p ; i < nsplit ; ++i ) {
715 Xsplit[i] = new double [_n];
716 for ( j = 0 ; j < _n ; ++j )
717 Xsplit[i][j] = _XX[i-_p][j];
718 }
719
720 // get the rectangle:
721 if ( _tgp_rect )
722 delete_matrix ( _tgp_rect );
723 _tgp_rect = get_data_rect ( Xsplit , nsplit , _n );
724
725 // TGP parameters:
726 Params tgp_params ( _n );
727 double * dparams = NOMAD::TGP_Model::get_TGP_dparams ( _n );
728 tgp_params.read_double ( dparams );
729 delete [] dparams;
730
731 int BTE[3];
732 if ( !get_BTE ( BTE ) ) {
733 for ( i = 0 ; i < nsplit ; ++i )
734 delete [] Xsplit[i];
735 delete [] Xsplit;
736 return false;
737 }
738
739 // display BTE:
740 #ifdef TGP_DEBUG
741 _out << std::endl
742 << "BTE={" << BTE[0] << ", " << BTE[1] << ", " << BTE[2] << "}"
743 << std::endl;
744 #endif
745
746 // compute the individual TGP models (one for each output):
747 double * ZZ = new double [_n_XX];
748
749 // Ds2x, expected reduction in predictive variance:
750 if ( _Ds2x == NULL ) {
751 _Ds2x = new double * [m];
752 for ( i = 0 ; i < m ; ++i )
753 _Ds2x[i] = NULL;
754 }
755 else {
756 for ( i = 0 ; i < m ; ++i )
757 if ( _Ds2x[i] ) {
758 delete [] _Ds2x[i];
759 _Ds2x[i] = NULL;
760 }
761 }
762
763 // improv, expected improvement of the objective (ranks):
764 if ( _improv ) {
765 delete [] _improv;
766 _improv = NULL;
767 }
768 if ( compute_improv )
769 _improv = new int [_n_XX];
770
771 for ( i = 0 ; i < m ; ++i ) {
772
773 if ( _tgp_models[i] ) {
774
775 _Ds2x[i] = ( compute_Ds2x ) ? new double [_n_XX] : NULL;
776
777 _tgp_models[i]->compute ( _X ,
778 _XX ,
779 Xsplit ,
780 _n ,
781 _n_XX ,
782 nsplit ,
783 &tgp_params ,
784 _tgp_rect ,
785 BTE ,
786 _tgp_linburn ,
787 _tgp_verb ,
788 ZZ ,
789 _Ds2x[i] ,
790 (i==index_obj) ? _improv : NULL );
791
792 // set XX_pts outputs #i:
793 for ( j = 0 ; j < _n_XX ; ++j ) {
794 XX_filtered_pts[j]->set_bb_output ( i , ZZ[j] );
795 XX_filtered_pts[j]->set_eval_status ( NOMAD::EVAL_OK );
796 }
797
798 // check if TGP has been interrupted:
799 if ( NOMAD::TGP_Output_Model::get_force_quit() ) {
800 _error_str = "TGP has been interrupted with ctrl-c";
801 break;
802 }
803 }
804 }
805
806 // clear memory:
807 for ( i = 0 ; i < nsplit ; ++i )
808 delete [] Xsplit[i];
809 delete [] Xsplit;
810 delete [] ZZ;
811
812 _model_computed = _error_str.empty();
813
814 return _model_computed;
815 }
816
817 /*--------------------------------------------*/
818 /* prediction at one point */
819 /* (x can be of size _n or _n0) */
820 /*--------------------------------------------*/
predict(NOMAD::Eval_Point & x,bool pred_outside_bnds)821 bool NOMAD::TGP_Model::predict ( NOMAD::Eval_Point & x , bool pred_outside_bnds )
822 {
823 if ( !_error_str.empty() )
824 return false;
825
826 if ( !_model_computed ) {
827 _error_str = "NOMAD::TGP_Model::compute() has not been called";
828 return false;
829 }
830
831 int i , i0 , ix , m = x.get_m() , nx = x.size();
832
833 // reset point outputs:
834 x.set_eval_status ( NOMAD::EVAL_FAIL );
835 for ( i = 0 ; i < m ; ++i )
836 x.set_bb_output ( i , NOMAD::Double() );
837
838 // check dimensions:
839 if ( m != static_cast<int>(_bbot.size()) ||
840 ( nx != _n0 && nx != _n ) ) {
841 _error_str = "predict error: bad x dimensions";
842 return false;
843 }
844
845 double ZZ , * XX = new double[_n];
846
847 // set the coordinates and check the bounds:
848 for ( i = 0 ; i < _n ; ++i ) {
849
850 ix = ( nx == _n0 ) ? _fv_index[i] : i;
851
852 if ( !pred_outside_bnds ) {
853 i0 = _fv_index[i];
854 if ( x[ix] < _lb[i0] || x[ix] > _ub[i0] ) {
855 delete [] XX;
856 return false; // this is not an error
857 }
858 }
859
860 XX[i] = x[ix].value();
861 }
862
863 // predictions (one for each output):
864 for ( i = 0 ; i < m ; ++i )
865 if ( _tgp_models[i] ) {
866 if ( !_tgp_models[i]->predict ( XX ,
867 _n ,
868 ZZ ,
869 _tgp_rect ) ) {
870 std::ostringstream oss;
871 oss << "predict error: problem with model #" << i;
872 _error_str = oss.str();
873 break;
874 }
875 x.set_bb_output ( i , ZZ );
876 }
877
878 delete [] XX;
879
880 if ( !_error_str.empty() ) {
881 x.set_eval_status ( NOMAD::EVAL_FAIL );
882 return false;
883 }
884
885 x.set_eval_status ( NOMAD::EVAL_OK );
886 return true;
887 }
888
889 /*-----------------------------------------------------------------*/
890 /* this function checks if the p x n matrix X is of full rank by */
891 /* applying the Cholesky decomposition to the sym. def. pos. */
892 /* nxn matrix X'X */
893 /* (static, private) */
894 /*-----------------------------------------------------------------*/
check_full_rank(double ** X,int p,int n)895 bool NOMAD::TGP_Model::check_full_rank ( double ** X , int p , int n ) {
896
897 int i , j , k , ki , kii , kij , kj , kji , nn12 = n*(n+1)/2;
898
899 // create XTX (X'X):
900 double * XTX = new double [nn12];
901 for ( i = 0 ; i < n ; ++i ) {
902 ki = i*(i+1)/2;
903 for ( j = 0 ; j <= i ; ++j ) {
904 kij = ki + j;
905 XTX[kij] = 0.0;
906 for ( k = 0 ; k < p ; ++k )
907 XTX[kij] += X[k][i] * X[k][j];
908 }
909 }
910
911 // create chol:
912 double * chol = new double [nn12] , tmp1 , tmp2 , eps = 1e-10;
913
914 // Choleski decomposition:
915 for ( i = 0 ; i < n ; ++i ) {
916
917 ki = i*(i+1)/2;
918 kii = ki+i;
919 tmp1 = XTX[kii];
920 for ( k = 0 ; k < i ; ++k )
921 tmp1 -= pow(chol[ki+k],2.0);
922 if ( fabs ( tmp1 ) <= eps ) {
923 delete [] XTX;
924 delete [] chol;
925 return false;
926 }
927
928 if ( i == n-1 )
929 break;
930
931 tmp1 = sqrt(tmp1);
932 chol[kii] = tmp1;
933
934 for ( j = i+1 ; j < n ; ++j ) {
935 kj = j*(j+1)/2;
936 kji = kj+i;
937 tmp2 = XTX[kji];
938 for ( k = 0 ; k < i ; ++k )
939 tmp2 -= chol[ki+k]*chol[kj+k];
940 chol[kji] = tmp2/tmp1;
941 }
942 }
943
944 delete [] XTX;
945 delete [] chol;
946
947 return true;
948 }
949
950 /*----------------------------------------------------------------*/
951 /* compute model h and f values given one blackbox output */
952 /*----------------------------------------------------------------*/
eval_hf(const NOMAD::Point & bbo,const NOMAD::Double & h_min,NOMAD::hnorm_type h_norm,NOMAD::Double & h,NOMAD::Double & f) const953 void NOMAD::TGP_Model::eval_hf ( const NOMAD::Point & bbo ,
954 const NOMAD::Double & h_min ,
955 NOMAD::hnorm_type h_norm ,
956 NOMAD::Double & h ,
957 NOMAD::Double & f ) const
958 {
959 f.clear();
960 h = 0.0;
961
962 int m = bbo.size();
963
964 if ( m != static_cast<int>(_bbot.size()) )
965 throw NOMAD::Exception ( "TGP_Model.cpp" , __LINE__ ,
966 "TGP_Model::eval_hf() called with an invalid bbo argument" );
967
968 NOMAD::Double bboi;
969
970 for ( int i = 0 ; i < m ; ++i ) {
971
972 bboi = bbo[i];
973
974 if ( bboi.is_defined() ) {
975
976 if ( _bbot[i] == NOMAD::EB || _bbot[i] == NOMAD::PEB_E ) {
977
978 if ( bboi > h_min ) {
979 h.clear();
980 return;
981 }
982 }
983
984 else if ( ( _bbot[i] == NOMAD::FILTER ||
985 _bbot[i] == NOMAD::PB ||
986 _bbot[i] == NOMAD::PEB_P ) ) {
987 if ( bboi > h_min ) {
988 switch ( h_norm ) {
989 case NOMAD::L1:
990 h += bboi;
991 break;
992 case NOMAD::L2:
993 h += bboi * bboi;
994 break;
995 case NOMAD::LINF:
996 if ( bboi > h )
997 h = bboi;
998 break;
999 }
1000 }
1001 }
1002
1003 else if ( _bbot[i] == NOMAD::OBJ )
1004 f = bboi;
1005 }
1006
1007 }
1008
1009 if ( h_norm == NOMAD::L2 )
1010 h = h.sqrt();
1011 }
1012
1013 /*---------------------------------------*/
1014 /* compute the TGP dparam array */
1015 /* (static, private) */
1016 /*---------------------------------------*/
get_TGP_dparams(int n)1017 double * NOMAD::TGP_Model::get_TGP_dparams ( int n ) {
1018
1019 double * dparams = new double [n*(n+3)+41];
1020 int i , j , k;
1021
1022 // tree (p <- c(as.numeric(params$tree))):
1023 dparams[0] = 0.5;
1024 dparams[1] = 2.0;
1025 dparams[2] = n+2;
1026 if ( dparams[2] < 10 )
1027 dparams[2] = 10;
1028 dparams[3] = 1.0;
1029 dparams[4] = n;
1030
1031 // params$meanfn == "linear" :
1032 dparams[5] = 0.0;
1033
1034 // params$bprior == "bflat" :
1035 dparams[6] = 2.0;
1036
1037 // p <- c(p, as.numeric(params$beta)) : n+1 zeros
1038 k = 6;
1039 for ( i = 0 ; i <= n ; ++i )
1040 dparams[++k] = 0.0;
1041
1042 // p <- c(p, as.numeric(params$Wi)) : I_{n+1}
1043 for ( i = 0 ; i <= n ; ++i )
1044 for ( j = 0 ; j <= n ; ++j )
1045 dparams[++k] = (i!=j) ? 0.0 : 1.0;
1046
1047 // p <- c(p, as.numeric(params$s2tau2)):
1048 dparams[++k] = 1.0;
1049 dparams[++k] = 1.0;
1050
1051 // p <- c(p, as.numeric(params$s2.p)):
1052 dparams[++k] = 5.0;
1053 dparams[++k] = 10.0;
1054
1055 // p <- c(p, as.numeric(params$s2.lam)):
1056 dparams[++k] = 0.2;
1057 dparams[++k] = 10.0;
1058
1059 // p <- c(p, as.numeric(params$tau2.p)):
1060 dparams[++k] = 5.0;
1061 dparams[++k] = 10.0;
1062
1063 // p <- c(p, as.numeric(params$tau2.lam)):
1064 dparams[++k] = 0.2;
1065 dparams[++k] = 0.1;
1066
1067 // params$corr == "expsep" :
1068 dparams[++k] = 1.0;
1069
1070 // p <- c(p, as.numeric(params$gd)):
1071 dparams[++k] = 0.1;
1072 dparams[++k] = 0.5;
1073
1074 // p <- c(p, as.numeric(params$nug.p)):
1075 dparams[++k] = 1.0;
1076 dparams[++k] = 1.0;
1077 dparams[++k] = 1.0;
1078 dparams[++k] = 1.0;
1079
1080 // if (params$nug.lam[1] == "fixed"), p <- c(p, rep(-1, 4)) :
1081 dparams[++k] = -1.0;
1082 dparams[++k] = -1.0;
1083 dparams[++k] = -1.0;
1084 dparams[++k] = -1.0;
1085
1086 // p <- c(p, as.numeric(params$gamma)):
1087 dparams[++k] = 0.0;
1088 dparams[++k] = 0.2;
1089 dparams[++k] = 0.7;
1090
1091 // p <- c(p, as.numeric(params$d.p)):
1092 dparams[++k] = 1.0;
1093 dparams[++k] = 20.0;
1094 dparams[++k] = 10.0;
1095 dparams[++k] = 10.0;
1096
1097 // if (params$d.lam[1] == "fixed"), p <- c(p, rep(-1, 4)):
1098 dparams[++k] = -1.0;
1099 dparams[++k] = -1.0;
1100 dparams[++k] = -1.0;
1101 dparams[++k] = -1.0;
1102
1103 return dparams;
1104 }
1105
1106 /*--------------------------------------------*/
1107 /* display the interpolation set X */
1108 /*--------------------------------------------*/
display_X(const NOMAD::Display & out,int display_limit) const1109 void NOMAD::TGP_Model::display_X ( const NOMAD::Display & out ,
1110 int display_limit ) const
1111 {
1112 if ( _p == 0 || !_X ) {
1113 out << "no interpolation points" << std::endl;
1114 return;
1115 }
1116
1117 int i , j;
1118 int m = _bbot.size();
1119 int i0 = ( display_limit > 0 ) ? _p - display_limit : 0;
1120 NOMAD::Point x(_n) , bbo(m);
1121
1122 out << NOMAD::open_block ( "interpolation points (X)");
1123
1124 if ( i0 > 0 )
1125 out << "..." << std::endl;
1126 else if ( i0 < 0 )
1127 i0 = 0;
1128
1129 for ( i = i0 ; i < _p ; ++i ) {
1130
1131 for ( j = 0 ; j < _n ; ++j )
1132 x[j] = _X[i][j];
1133
1134 bbo = NOMAD::Point(m);
1135 if ( _tgp_models )
1136 for ( j = 0 ; j < m ; ++j )
1137 if ( _tgp_models[j] )
1138 bbo[j] = (_tgp_models[j]->get_Z())[i];
1139
1140 out << "#";
1141 out.display_int_w ( i , _p );
1142 out << " x=(";
1143 x.display ( out , " " , 15 , -1 );
1144 out << " ) f(x)=[";
1145 bbo.display ( out , " " , 15 , -1 );
1146 out << " ]" << std::endl;
1147 }
1148
1149 std::ostringstream oss;
1150 oss << "(size=" << _p << ")";
1151 out << NOMAD::close_block ( oss.str() ) << std::endl;
1152 }
1153
1154 /*---------------------------------------------------------*/
1155 /* get the XX points with the largest expected reduction */
1156 /* in predictive variance, for each output */
1157 /* (no duplicates) */
1158 /*---------------------------------------------------------*/
get_Ds2x_points(std::set<int> & pts_indexes) const1159 void NOMAD::TGP_Model::get_Ds2x_points ( std::set<int> & pts_indexes ) const
1160 {
1161 pts_indexes.clear();
1162 if ( !_Ds2x || _n_XX == 0 )
1163 return;
1164
1165 int i , j , k , m = _bbot.size();
1166
1167 for ( i = 0 ; i < m ; ++i )
1168 if ( _Ds2x[i] ) {
1169
1170 NOMAD::Double max;
1171 k = -1;
1172 for ( j = 0 ; j < _n_XX ; ++j )
1173 if ( !max.is_defined() || _Ds2x[i][j] > max ) {
1174 max = _Ds2x[i][j];
1175 k = j;
1176 }
1177
1178 if ( k >= 0 )
1179 pts_indexes.insert ( k );
1180 }
1181 }
1182
1183 /*----------------------------------------------------------------------------*/
1184 /* get the XX points with the largest expected improvement of the objective */
1185 /*----------------------------------------------------------------------------*/
get_improv_points(std::list<int> & pts_indexes) const1186 void NOMAD::TGP_Model::get_improv_points ( std::list<int> & pts_indexes ) const
1187 {
1188 pts_indexes.clear();
1189 if ( !_improv || _n_XX == 0 )
1190 return;
1191
1192 int j;
1193 NOMAD::Point * XX_pt;
1194 std::multiset<NOMAD::Model_Sorted_Point> pts;
1195 std::multiset<NOMAD::Model_Sorted_Point>::const_iterator it , end;
1196
1197 // 1. sort:
1198 for ( j = 0 ; j < _n_XX ; ++j ) {
1199 XX_pt = new NOMAD::Point ( 1 );
1200 (*XX_pt)[0] = j;
1201 pts.insert ( NOMAD::Model_Sorted_Point ( XX_pt , _improv[j] ) );
1202 }
1203
1204 // 2. construct pts_indexes (exclude points with improv >= n_XX):
1205 end = pts.end();
1206 for ( it = pts.begin() ; it != end ; ++it ) {
1207 if ( it->get_dist() < _n_XX )
1208 pts_indexes.push_back ( static_cast<int> ( (*it->get_point())[0].value() ) );
1209 delete it->get_point();
1210 }
1211 }
1212
1213 /*--------------------------------------------------------------*/
1214 /* display the expected improvement of the objective (improv) */
1215 /* (better ranks are displayed first) */
1216 /*--------------------------------------------------------------*/
display_improv(const NOMAD::Display & out,int display_limit) const1217 void NOMAD::TGP_Model::display_improv ( const NOMAD::Display & out ,
1218 int display_limit ) const
1219 {
1220 if ( !_improv || _n_XX == 0 ) {
1221 out << "improv has not been computed" << std::endl;
1222 return;
1223 }
1224
1225 int j , k;
1226 NOMAD::Point * XX_pt;
1227 std::multiset<NOMAD::Model_Sorted_Point> pts;
1228 std::multiset<NOMAD::Model_Sorted_Point>::const_iterator it , end;
1229
1230 // 1. sort:
1231 for ( j = 0 ; j < _n_XX ; ++j ) {
1232
1233 // construct a NOMAD::Point from matrix _XX:
1234 XX_pt = new NOMAD::Point ( _fv );
1235 for ( k = 0 ; k < _n0 ; ++k )
1236 if ( _av_index[k] >= 0 )
1237 (*XX_pt)[k] = _XX[j][_av_index[k]];
1238
1239 // insert this point in the sorted list:
1240 pts.insert ( NOMAD::Model_Sorted_Point ( XX_pt , _improv[j] ) );
1241 }
1242
1243 // 2. display:
1244 end = pts.end();
1245 for ( j = 0 , it = pts.begin() ; it != end ; ++it , ++j ) {
1246
1247 if ( display_limit <= 0 || j < display_limit ) {
1248 out << "x=( ";
1249 it->get_point()->display ( out , " " , 6 , -1 );
1250 out << " ) improv=" << it->get_dist() << std::endl;
1251 }
1252
1253 else if ( j == display_limit )
1254 out << "..." << std::endl;
1255
1256 delete it->get_point();
1257 }
1258 }
1259
1260 /*----------------------------------------------------------------*/
1261 /* display the expected reduction in predictive variance (Ds2x) */
1262 /* (larger values are displayed first) */
1263 /*----------------------------------------------------------------*/
display_Ds2x(const NOMAD::Display & out,int display_limit) const1264 void NOMAD::TGP_Model::display_Ds2x ( const NOMAD::Display & out ,
1265 int display_limit ) const
1266 {
1267 if ( !_Ds2x || _n_XX == 0 ) {
1268 out << "matrix Ds2x has not been computed" << std::endl;
1269 return;
1270 }
1271
1272 int i , j , k , m = _bbot.size();
1273 NOMAD::Point * XX_pt;
1274 std::multiset<NOMAD::Model_Sorted_Point> pts;
1275 std::multiset<NOMAD::Model_Sorted_Point>::const_iterator it , end;
1276
1277 for ( i = 0 ; i < m ; ++i ) {
1278
1279 if ( m > 1 ) {
1280 std::ostringstream oss;
1281 oss << "output #" << i;
1282 out << NOMAD::open_block ( oss.str() );
1283 }
1284
1285 if ( _Ds2x[i] ) {
1286
1287 // 1. sort:
1288 for ( j = 0 ; j < _n_XX ; ++j ) {
1289
1290 // construct a NOMAD::Point from matrix _XX:
1291 XX_pt = new NOMAD::Point ( _fv );
1292 for ( k = 0 ; k < _n0 ; ++k )
1293 if ( _av_index[k] >= 0 )
1294 (*XX_pt)[k] = _XX[j][_av_index[k]];
1295
1296 // insert this point in the sorted list:
1297 pts.insert ( NOMAD::Model_Sorted_Point ( XX_pt , -_Ds2x[i][j] ) );
1298 }
1299
1300 // 2. display:
1301 end = pts.end();
1302 for ( j = 0 , it = pts.begin() ; it != end ; ++it , ++j ) {
1303
1304 if ( display_limit <= 0 || j < display_limit ) {
1305 out << "x=( ";
1306 it->get_point()->display ( out , " " , 6 , -1 );
1307 out << " ) Ds2x=" << it->get_dist()*-1.0 << std::endl;
1308 }
1309
1310 else if ( j == display_limit )
1311 out << "..." << std::endl;
1312
1313 delete it->get_point();
1314 }
1315
1316 pts.clear();
1317 }
1318 else
1319 out << "NULL" << std::endl;
1320
1321 if ( m > 1 )
1322 out.close_block();
1323 }
1324 }
1325
1326 /*-------------------------------------------------------*/
1327 /* display the error stats for the interpolation set X */
1328 /*-------------------------------------------------------*/
display_X_errors(const NOMAD::Display & out)1329 void NOMAD::TGP_Model::display_X_errors ( const NOMAD::Display & out )
1330 {
1331 if ( _p == 0 || !_X ) {
1332 out << "no interpolation points" << std::endl;
1333 return;
1334 }
1335
1336 int i , j , m = _bbot.size();
1337 NOMAD::Point min_err(m) , max_err(m) , avg_err(m,0.0) , sd_err(m,0.0);
1338 NOMAD::Eval_Point x ( _n , m );
1339 double ** err = new double * [_p];
1340
1341 for ( i = 0 ; i < _p ; ++i ) {
1342
1343 err[i] = new double[m];
1344
1345 for ( j = 0 ; j < _n ; ++j )
1346 x[j] = _X[i][j];
1347
1348 if ( predict ( x , true ) ) {
1349
1350 for ( j = 0 ; j < m ; ++j )
1351 if ( _tgp_models[j] ) {
1352
1353 // relative error (in %) for point #i and output #j:
1354 err[i][j] = ( x.get_bb_outputs()[j].rel_err((_tgp_models[j]->get_Z())[i])
1355 * 100.0).value();
1356
1357 // out << "f=" << (_tgp_models[j]->get_Z())[i] << " "
1358 // << "m=" << x.get_bb_outputs()[j] << " err=" << err[i][j]
1359 // << std::endl;
1360
1361 if ( !min_err[j].is_defined() || err[i][j] < min_err[j].value() )
1362 min_err[j] = err[i][j];
1363
1364 if ( !max_err[j].is_defined() || err[i][j] > max_err[j].value() )
1365 max_err[j] = err[i][j];
1366
1367 avg_err[j] += err[i][j];
1368 }
1369 }
1370 else {
1371 for ( j = 0 ; j <= i ; ++j )
1372 delete [] err[j];
1373 delete [] err;
1374 out << "cannot predict interpolation errors ("
1375 << _error_str << ")" << std::endl;
1376 return;
1377 }
1378 }
1379
1380 for ( j = 0 ; j < m ; ++j )
1381 if ( _tgp_models[j] ) {
1382
1383 // compute the median error:
1384 NOMAD::Double med_err;
1385 {
1386 if ( _p == 1 )
1387 med_err = err[0][j];
1388 else if ( _p == 2 )
1389 med_err = ( err[0][j] + err[1][j] ) / 2.0;
1390
1391 else {
1392 std::multiset<double> sorted_errors;
1393 for ( i = 0 ; i < _p ; ++i )
1394 sorted_errors.insert ( err[i][j] );
1395 std::multiset<double>::const_iterator it , end = sorted_errors.end();
1396 --end;
1397 for ( it = sorted_errors.begin() , i = 0 ; it != end ; ++it , ++i ) {
1398 if ( i == (_p+1)/2-1 ) {
1399 med_err = *it;
1400 if ( _p%2==0 ) {
1401 ++it;
1402 med_err = ( med_err + *it ) / 2.0;
1403 }
1404 break;
1405 }
1406 }
1407 }
1408 }
1409
1410 // compute the mean and the standard deviation:
1411 avg_err[j] /= _p;
1412 for ( i = 0 ; i < _p ; ++i )
1413 sd_err[j] += ( avg_err[j] - err[i][j] ).pow2();
1414 sd_err[j] = (sd_err[j] / _p).sqrt();
1415
1416 // display:
1417 if ( m > 1 ) {
1418 std::ostringstream oss;
1419 oss << "output #" << j;
1420 if ( _tgp_models[j]->is_fixed() )
1421 oss << " (fixed)";
1422 else if ( _tgp_models[j]->is_binary() )
1423 oss << " (binary)";
1424 out << NOMAD::open_block ( oss.str() );
1425 }
1426
1427 out << "min : ";
1428 min_err[j].display ( out , "%6.2f" );
1429 out << std::endl << "max : ";
1430 max_err[j].display ( out , "%6.2f" );
1431 out << std::endl << "median: ";
1432 med_err.display ( out , "%6.2f" );
1433 out << std::endl << "mean : ";
1434 avg_err[j].display ( out , "%6.2f" );
1435 out << std::endl << "sd : ";
1436 sd_err[j].display ( out , "%6.2f" );
1437 out << std::endl;
1438
1439 if ( m > 1 )
1440 out.close_block();
1441 }
1442
1443 for ( i = 0 ; i < _p ; ++i )
1444 delete [] err[i];
1445 delete [] err;
1446 }
1447
1448 /*--------------------------------------------*/
1449 /* display */
1450 /*--------------------------------------------*/
display(const NOMAD::Display & out) const1451 void NOMAD::TGP_Model::display ( const NOMAD::Display & out ) const
1452 {
1453 if ( !_error_str.empty() ) {
1454 out << "error with model" << std::endl;
1455 return;
1456 }
1457
1458 int i , j;
1459
1460 // fixed variables:
1461 out << "fixed_var = [ " << _fv << "]" << std::endl;
1462 out << "av_index = ";
1463 if ( _av_index ) {
1464 out << "[ ";
1465 for ( i = 0 ; i < _n0 ; ++i )
1466 out << _av_index[i] << " ";
1467 out << "]" << std::endl;
1468 }
1469 else
1470 out << "NULL" << std::endl;
1471 out << "fv_index = ";
1472 if ( _fv_index ) {
1473 out << "[ ";
1474 for ( i = 0 ; i < _n ; ++i )
1475 out << _fv_index[i] << " ";
1476 out << "]" << std::endl;
1477 }
1478 else
1479 out << "NULL" << std::endl;
1480
1481 // bounds:
1482 out << "lb = [ " << _lb << "]" << std::endl
1483 << "ub = [ " << _ub << "]" << std::endl
1484 << std::endl;
1485
1486 // display X:
1487 if ( !_X )
1488 out << "X = NULL" << std::endl;
1489 else {
1490 out << "X = [";
1491 for ( i = 0 ; i < _p ; ++i ) {
1492 out << "\t";
1493 for ( j = 0 ; j < _n ; ++j )
1494 out << std::setw(15) << _X[i][j] << " ";
1495 out << ( (i==_p-1) ? "]" : ";" ) << std::endl;
1496 }
1497 out << "size(X)=" << _p << "x" << _n << std::endl << std::endl;
1498 }
1499
1500 // display XX:
1501 if ( !_XX )
1502 out << "XX = NULL" << std::endl;
1503 else {
1504 out << "XX = [";
1505 for ( i = 0 ; i < _n_XX ; ++i ) {
1506 out << "\t";
1507 for ( j = 0 ; j < _n ; ++j )
1508 out << std::setw(15) << _XX[i][j] << " ";
1509 out << ( (i==_n_XX-1) ? "]" : ";" ) << std::endl;
1510 }
1511 out << "size(XX)=" << _n_XX << "x" << _n << std::endl << std::endl;
1512 }
1513
1514 // display models:
1515 out << std::endl;
1516 if ( _tgp_models ) {
1517 int m = _bbot.size();
1518 for ( i = 0 ; i < m ; ++i ) {
1519 if ( _tgp_models[i] ) {
1520 std::ostringstream oss;
1521 oss << "model #" << i;
1522 out.open_block ( oss.str() );
1523 _tgp_models[i]->display ( out );
1524 out.close_block();
1525 }
1526 else
1527 out << "model #" << i << ": NULL" << std::endl;
1528 out << std::endl;
1529 }
1530 }
1531 else
1532 out << "no models" << std::endl << std::endl;
1533 }
1534
1535 /*----------------------------------------------------------------*/
1536 /* access to the width of the interpolation set */
1537 /*----------------------------------------------------------------*/
1538 #ifdef MODEL_STATS
get_Yw(void) const1539 NOMAD::Double NOMAD::TGP_Model::get_Yw ( void ) const
1540 {
1541 NOMAD::Double Yw , tmp;
1542 for ( int i = 0 ; i < _n0 ; ++i ) {
1543 tmp = _ub[i]-_lb[i];
1544 if ( !Yw.is_defined() || tmp > Yw )
1545 Yw = tmp;
1546 }
1547 return Yw;
1548 }
1549 #endif
1550
1551 #endif
1552