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_Output_Model.cpp
38 \brief TGP (Bayesian treed Gaussian process) model for one output (implementation)
39 \author Sebastien Le Digabel
40 \date 2011-02-07
41 \see TGP_Output_Model.hpp
42 */
43
44 #ifndef USE_TGP
45
46 int TGP_OUTPUT_MODEL_DUMMY; // avoids that TGP_Output_Model.o has no symbols with ranlib
47
48 #else
49
50 #include "TGP_Output_Model.hpp"
51
52 /*---------------------------------------------------------*/
53 /* NOMAD-TGP callback function (called regularly by TGP) */
54 /*---------------------------------------------------------*/
55 // SLD -- 2014-09-04
56 // void NOMAD::TGP_callback ( bool & TGP_interrupt )
57 // {
58 // if ( NOMAD::TGP_Output_Model::get_force_quit() )
59 // TGP_interrupt = true;
60 // }
61
62 /*-----------------------------------*/
63 /* static members initialization */
64 /*-----------------------------------*/
65 double NOMAD::TGP_Output_Model::_ditemps[7] =
66 { 1.0 , 0.0 , 0.0 , 1.0 , 1.0 , 0.0 , 1.0 };
67
68 bool NOMAD::TGP_Output_Model::_force_quit = false;
69
70 /*-----------------------------------*/
71 /* constructor */
72 /*-----------------------------------*/
TGP_Output_Model(const std::list<const NOMAD::Eval_Point * > & X_pts,int bbo_index,int seed,const NOMAD::Display & out)73 NOMAD::TGP_Output_Model::TGP_Output_Model
74 ( const std::list<const NOMAD::Eval_Point *> & X_pts ,
75 int bbo_index ,
76 int seed ,
77 const NOMAD::Display & out )
78 : _out ( out ) ,
79 _p ( X_pts.size() ) ,
80 _Z ( new double[_p] ) ,
81 _Z_is_scaled ( false ) ,
82 _is_binary ( true ) ,
83 _bin_values ( 2 ) ,
84 _is_fixed ( false ) ,
85 _tgp_state ( NULL ) ,
86 _tgp_model ( NULL ) ,
87 _tgp_its ( NULL )
88 {
89 NOMAD::TGP_Output_Model::_force_quit = false;
90
91 _Z_scaling[0] = _Z_scaling[1] = 0.0;
92
93 std::list<const NOMAD::Eval_Point *>::const_iterator it , end = X_pts.end();
94
95 NOMAD::Double tmp , Zmin , Zmax , Zsum = 0.0;
96 int j = 0;
97
98 for ( it = X_pts.begin() ; it != end ; ++it ) {
99
100 // the output value:
101 tmp = (*it)->get_bb_outputs()[bbo_index];
102 _Z[j++] = tmp.value();
103
104 // Z scaling parameters (1/2):
105 Zsum += tmp;
106 if ( !Zmin.is_defined() || tmp < Zmin )
107 Zmin = tmp;
108 if ( !Zmax.is_defined() || tmp > Zmax )
109 Zmax = tmp;
110
111 // check if the output is binary:
112 if ( _is_binary ) {
113 if ( !_bin_values[0].is_defined() )
114 _bin_values[0] = tmp;
115 else if ( !_bin_values[1].is_defined() && tmp != _bin_values[0] )
116 _bin_values[1] = tmp;
117 else {
118 if ( tmp != _bin_values[0] && tmp != _bin_values[1] )
119 _is_binary = false;
120 }
121 }
122 }
123
124 // Z scaling parameters (1/2):
125 _Z_scaling[0] = (Zmax - Zmin).value();
126
127 // the output is fixed:
128 if ( _Z_scaling[0] == 0.0 )
129 _is_fixed = true;
130
131 else {
132
133 _Z_scaling[1] = (Zsum/_p).value() / _Z_scaling[0];
134
135 if ( !_is_binary )
136 _bin_values = NOMAD::Point(2);
137
138 // RNG (random number generator):
139 int state[] = { 896 , 265 , 371 };
140
141 // if seed==0, the model will be the same as the one constructed in R,
142 // with values taken from the R tgp functions for:
143 // set.seed(0)
144 // state <- sample(seq(0, 999), 3)
145
146 // otherwise use rand() to get three different integers in [0;999]:
147 if ( seed != 0 ) {
148 state[0] = rand()%1000;
149 while ( state[1] == state[0] )
150 state[1] = rand()%1000;
151 while ( state[2] == state[0] || state[2] == state[1] )
152 state[2] = rand()%1000;
153 }
154 _tgp_state = newRNGstate ( three2lstate ( state ) );
155
156 // importance tempering:
157 _tgp_its = new Temper ( NOMAD::TGP_Output_Model::_ditemps );
158 }
159 }
160
161 /*--------------------------------------------*/
162 /* destructor */
163 /*--------------------------------------------*/
~TGP_Output_Model(void)164 NOMAD::TGP_Output_Model::~TGP_Output_Model ( void )
165 {
166 if ( _Z )
167 delete [] _Z;
168 if ( _tgp_model )
169 delete _tgp_model;
170 if ( _tgp_its )
171 delete _tgp_its;
172 if ( _tgp_state )
173 deleteRNGstate ( _tgp_state );
174 }
175
176 /*--------------------------------------------*/
177 /* compute the model */
178 /*--------------------------------------------*/
compute(double ** X,double ** XX,double ** Xsplit,int n,int n_XX,int nsplit,Params * tgp_params,double ** tgp_rect,int * tgp_BTE,bool tgp_linburn,bool tgp_verb,double * ZZ,double * Ds2x,int * improv)179 void NOMAD::TGP_Output_Model::compute ( double ** X ,
180 double ** XX ,
181 double ** Xsplit ,
182 int n ,
183 int n_XX ,
184 int nsplit ,
185 Params * tgp_params ,
186 double ** tgp_rect ,
187 int * tgp_BTE ,
188 bool tgp_linburn ,
189 bool tgp_verb ,
190 double * ZZ , // OUT
191 double * Ds2x , // OUT, may be NULL
192 int * improv ) // OUT, may be NULL
193 {
194 bool compute_Ds2x = ( Ds2x != NULL );
195 bool compute_improv = ( improv != NULL );
196
197 // the output is fixed:
198 if ( _is_fixed ) {
199 for ( int j = 0 ; j < n_XX ; ++j ) {
200 ZZ[j] = _bin_values[0].value();
201 if ( compute_Ds2x )
202 Ds2x[j] = 0.0;
203 }
204 return;
205 }
206
207 // scale Z:
208 scale_Z();
209
210 // construct the TGP model:
211 _tgp_model = new Model ( tgp_params ,
212 n ,
213 tgp_rect ,
214 0 , // Id=0
215 false , // trace=false
216 _tgp_state );
217
218 _tgp_model->Init ( X ,
219 _p ,
220 n ,
221 _Z ,
222 _tgp_its ,
223 NULL , // dtree=NULL
224 0 , // ncol=0
225 NULL ); // dhier=NULL
226
227 // set the NOMAD-TGP callback function:
228 // _tgp_model->set_callback ( NOMAD::TGP_callback ); // SLD -- 2014-09-04
229
230 // TGP verbim (display):
231 #ifdef TGP_DEBUG
232 _tgp_model->Outfile ( stdout , 1 ); // set 10 for full display
233 #else
234 if ( tgp_verb )
235 _tgp_model->Outfile ( stdout , 1 );
236 else
237 _tgp_model->Outfile ( NULL , 0 );
238 #endif
239
240 // set the splitting locations (Xsplit):
241 _tgp_model->set_Xsplit ( Xsplit , nsplit , n );
242
243 // linear model initialization rounds -B thru 1:
244 if ( tgp_linburn )
245 _tgp_model->Linburn ( tgp_BTE[0] , _tgp_state );
246
247 // do model rounds 1 thru B (burn in):
248 _tgp_model->Burnin ( tgp_BTE[0] , _tgp_state );
249
250 // do the MCMC rounds B,...,T:
251 Preds * tgp_preds = new_preds ( XX ,
252 n_XX ,
253 0 ,
254 n ,
255 tgp_rect ,
256 tgp_BTE[1]-tgp_BTE[0] ,
257 false , // pred_n
258 true , // krige
259 _tgp_its->IT_ST_or_IS() ,
260 compute_Ds2x ,
261 compute_improv ,
262 false , // sens
263 tgp_BTE[2] );
264
265 _tgp_model->Sample ( tgp_preds , tgp_BTE[1]-tgp_BTE[0] , _tgp_state );
266
267 // if importance tempering, then update the pseudo-prior
268 // based on the observation counts:
269 if ( _tgp_its->Numit() > 1 )
270 _tgp_its->UpdatePrior ( _tgp_model->update_tprobs() , _tgp_its->Numit() );
271
272 // copy back the itemps:
273 _tgp_model->DupItemps ( _tgp_its );
274
275 // kriging mean:
276 wmean_of_columns ( ZZ ,
277 tgp_preds->ZZm ,
278 tgp_preds->R ,
279 n_XX ,
280 (_tgp_its->IT_ST_or_IS()) ? tgp_preds->w : NULL );
281
282 int i;
283
284 // expected reduction in predictive variance (Ds2x):
285 if ( compute_Ds2x ) {
286 for ( i = 0 ; i < n_XX ; ++i )
287 Ds2x[i] = 0.0;
288 if ( tgp_preds->Ds2x )
289 wmean_of_columns ( Ds2x ,
290 tgp_preds->Ds2x ,
291 tgp_preds->R ,
292 tgp_preds->nn ,
293 (_tgp_its->IT_ST_or_IS()) ? tgp_preds->w : NULL );
294 }
295
296 // expected improvement of objective (improv):
297 if ( compute_improv ) {
298
299 // double * improvec = new double [n_XX];
300 // wmean_of_columns ( improvec ,
301 // tgp_preds->improv,
302 // tgp_preds->R,
303 // tgp_preds->nn,
304 // (_tgp_its->IT_ST_or_IS()) ? tgp_preds->w : NULL );
305 // for ( i = 0 ; i < n_XX ; ++i )
306 // _out << "IMPROVEC[" << i<< "] = " << improvec[i] << std::endl;
307 // delete [] improvec;
308
309 int *ir = (int*) GetImprovRank ( tgp_preds->R ,
310 tgp_preds->nn ,
311 tgp_preds->improv ,
312 true , // improv=true
313 tgp_preds->nn , // numirank = n_XX
314 (_tgp_its->IT_ST_or_IS()) ? tgp_preds->w : NULL );
315
316 for ( i = 0 ; i < n_XX ; ++i ) {
317 improv[i] = ir[i];
318 if ( improv[i] == 0 )
319 improv[i] = n_XX;
320 }
321
322 free ( ir );
323
324 // for ( i = 0 ; i < n_XX ; ++i )
325 // _out << "RANK[" << i<< "] = " << improv[i] << std::endl;
326 }
327
328 delete_preds ( tgp_preds );
329
330 #ifdef TGP_DEBUG
331 _tgp_model->Outfile ( NULL , 0 );
332 #endif
333
334 // unscale Z and ZZ:
335 unscale_Z ( ZZ , n_XX );
336 unscale_Z();
337
338 // treat binary output:
339 if ( _is_binary )
340 treat_binary_output ( ZZ , n_XX );
341
342 // disable TGP display:
343 _tgp_model->Outfile ( NULL , 0 );
344 }
345
346 /*--------------------------------------------*/
347 /* prediction at one point */
348 /*--------------------------------------------*/
predict(double * XX,int n,double & ZZ,double ** tgp_rect) const349 bool NOMAD::TGP_Output_Model::predict ( double * XX ,
350 int n ,
351 double & ZZ ,
352 double ** tgp_rect ) const
353 {
354 if ( _is_fixed ) {
355 ZZ = _bin_values[0].value();
356 return true;
357 }
358
359 // do the MCMC rounds B,...,T:
360 Preds * tgp_preds = new_preds ( &XX ,
361 1 ,
362 0 ,
363 n ,
364 tgp_rect ,
365 1 , // instead of T-B
366 false , // pred_n
367 true , // krige
368 _tgp_its->IT_ST_or_IS() ,
369 false , // delta_s2 (flag for Ds2x)
370 false , // improv
371 false , // sens
372 1 ); // instead of E
373
374 // new TGP function for the one point prediction:
375 _tgp_model->MAPreplace();
376
377 // prediction:
378 _tgp_model->Predict ( tgp_preds ,
379 1 , // instead of T-B
380 _tgp_state );
381
382 // kriging mean:
383 ZZ = *tgp_preds->ZZm[0];
384 // no need to do:
385 // wmean_of_columns ( &ZZ ,
386 // tgp_preds->ZZm ,
387 // tgp_preds->R ,
388 // 1 ,
389 // (_tgp_its->IT_ST_or_IS()) ? tgp_preds->w : NULL );
390
391 delete_preds ( tgp_preds );
392
393 // unscale:
394 unscale_Z ( &ZZ , 1 );
395
396 // treat binary output:
397 if ( _is_binary )
398 treat_binary_output ( &ZZ , 1 );
399
400 return true;
401 }
402
403 /*--------------------------------------------*/
404 /* scale Z or ZZ (private) */
405 /*--------------------------------------------*/
scale_Z(void)406 void NOMAD::TGP_Output_Model::scale_Z ( void )
407 {
408 if ( _Z_is_scaled || _is_fixed )
409 return;
410 scale_Z ( _Z , _p );
411 _Z_is_scaled = true;
412 }
413
scale_Z(double * Z,int n) const414 void NOMAD::TGP_Output_Model::scale_Z ( double * Z , int n ) const
415 {
416 if ( _is_fixed )
417 return;
418 for ( int i = 0 ; i < n ; ++i )
419 Z[i] = Z[i] / _Z_scaling[0] - _Z_scaling[1];
420 }
421
422 /*--------------------------------------------*/
423 /* unscale Z or ZZ (private) */
424 /*--------------------------------------------*/
unscale_Z(void)425 void NOMAD::TGP_Output_Model::unscale_Z ( void )
426 {
427 if ( !_Z_is_scaled || _is_fixed )
428 return;
429 unscale_Z ( _Z , _p );
430 _Z_is_scaled = false;
431 }
432
unscale_Z(double * Z,int n) const433 void NOMAD::TGP_Output_Model::unscale_Z ( double * Z , int n ) const
434 {
435 if ( _is_fixed )
436 return;
437 for ( int i = 0 ; i < n ; ++i )
438 Z[i] = ( Z[i] + _Z_scaling[1] ) * _Z_scaling[0];
439 }
440
441 /*--------------------------------------------*/
442 /* treat binary output (private) */
443 /*--------------------------------------------*/
treat_binary_output(double * ZZ,int nout) const444 void NOMAD::TGP_Output_Model::treat_binary_output ( double * ZZ , int nout ) const
445 {
446 NOMAD::Double d0 , d1;
447 for ( int j = 0 ; j < nout ; ++j ) {
448 d0 = (_bin_values[0] - ZZ[j]).abs();
449 d1 = (_bin_values[1] - ZZ[j]).abs();
450 if ( d0 < d1 )
451 ZZ[j] = _bin_values[0].value();
452 else
453 ZZ[j] = _bin_values[1].value();
454 }
455 }
456
457 /*--------------------------------------------*/
458 /* display */
459 /*--------------------------------------------*/
display(const NOMAD::Display & out) const460 void NOMAD::TGP_Output_Model::display ( const NOMAD::Display & out ) const
461 {
462 out << "Z (";
463 if ( !_Z_is_scaled )
464 out << "un";
465 out << "scaled)" << " = [ ";
466 for ( int i = 0 ; i < _p ; ++i ) {
467 out << std::setw(15) << _Z[i];
468 if ( (i+1)%5 == 0 )
469 out << " ;" << std::endl << " ";
470 else
471 out << " ";
472 }
473 out << "]" << std::endl
474 << "size(Z)=" << _p << std::endl;
475 if ( _is_fixed )
476 out << "fixed output";
477 else if ( _is_binary )
478 out << "binary output: values in {"
479 << _bin_values[0] << "," << _bin_values[1]
480 << "}";
481 out << std::endl
482 << "Z_scaling=[" << _Z_scaling[0] << "," << _Z_scaling[1]
483 << "]" << std::endl;
484 }
485
486 #endif
487