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 Quad_Model.cpp
38 \brief Quadratic regression or MFN interpolation model (implementation)
39 \author Sebastien Le Digabel
40 \date 2010-08-31
41 \see Quad_Model.hpp
42 */
43 #include "Quad_Model.hpp"
44
45 /*-----------------------------------------------------------*/
46 /* constructor */
47 /*-----------------------------------------------------------*/
Quad_Model(const NOMAD::Display & out,const std::vector<NOMAD::bb_output_type> & bbot,const NOMAD::Cache & cache,const NOMAD::Signature & signature)48 NOMAD::Quad_Model::Quad_Model
49 ( const NOMAD::Display & out ,
50 const std::vector<NOMAD::bb_output_type> & bbot ,
51 const NOMAD::Cache & cache ,
52 const NOMAD::Signature & signature )
53 : _out ( out ) ,
54 _bbot ( bbot ) ,
55 _interpolation_type ( NOMAD::UNDEFINED_INTERPOLATION_TYPE ) ,
56 _n ( signature.get_n() ) ,
57 _nfree ( _n ) ,
58 _fixed_vars ( new bool [_n] ) ,
59 _index ( NULL ) ,
60 _alpha ( NULL ) ,
61 _cache ( cache ) ,
62 _signature ( signature ) ,
63 _error_flag ( true )
64 {
65 for ( int i = 0 ; i < _n ; ++i )
66 _fixed_vars[i] = false;
67 init_alpha();
68 }
69
70 /*-----------------------------------------------------------*/
71 /* destructor */
72 /*-----------------------------------------------------------*/
~Quad_Model(void)73 NOMAD::Quad_Model::~Quad_Model ( void )
74 {
75 int m = static_cast<int> ( _bbot.size() );
76 for ( int i = 0 ; i < m ; ++i )
77 delete _alpha[i];
78 delete [] _alpha;
79 delete [] _fixed_vars;
80 delete [] _index;
81
82 // clear the interpolation set Y:
83 for ( size_t k = 0 ; k < _Y.size() ; ++k )
84 delete _Y[k];
85 }
86
87 /*-----------------------------------------------------------*/
88 /* initialize alpha, the model parameters (private) */
89 /*-----------------------------------------------------------*/
init_alpha(void)90 void NOMAD::Quad_Model::init_alpha ( void )
91 {
92 _n_alpha = ( _nfree + 1 ) * ( _nfree + 2 ) / 2;
93
94 int i , m = static_cast<int> ( _bbot.size() );
95
96 // initialize _alpha:
97 // ------------------
98 if ( _alpha )
99 {
100 for ( i = 0 ; i < m ; ++i )
101 delete _alpha[i];
102 delete [] _alpha;
103 }
104
105 _alpha = new NOMAD::Point * [m];
106
107 for ( i = 0 ; i < m ; ++i )
108 _alpha[i] = ( _bbot[i] == NOMAD::OBJ || NOMAD::bbot_is_constraint(_bbot[i]) ) ?
109 new NOMAD::Point ( _n_alpha ) : NULL;
110
111 // initialize _index:
112 // ------------------
113
114 // example: with 3 variables (X,Y,Z) with Y fixed.
115 // --------
116 // the problem is reduced to the two variables x=X and y=Z,
117 // and _index is corresponds to:
118 //
119 // 0 1 0 1 : index[0] = 0
120 // 1 X 1 x : index[1] = 1
121 // 2 Y 2 y : index[2] = 3
122 // 3 Z 3 .5 x^2 : index[3] = 4
123 // 4 .5 X^2 4 .5 y^2 : index[4] = 6
124 // 5 .5 Y^2 5 xy : index[5] = 8
125 // 6 .5 Z^2
126 // 7 XY
127 // 8 XZ
128 // 9 YZ
129 //
130 // If there are no fixed variables, index is of size (n+1)(n+2)/2
131 // with index[i] = i.
132
133 delete [] _index;
134
135 _index = new int [_n_alpha];
136
137 int nm1 = _n - 1;
138 int c1 = 2*_n + 1;
139 int c2 = 1;
140 int k1 , k2;
141
142 _index[0] = 0;
143 for ( i = 0 ; i < _n ; ++i )
144 {
145 if ( !_fixed_vars[i] )
146 {
147 _index[c2 ] = i+1;
148 _index[c2+_nfree] = i+1+_n;
149 ++c2;
150 }
151 }
152
153 c2 += _nfree;
154
155 for ( k1 = 0 ; k1 < nm1 ; ++k1 )
156 for ( k2 = k1+1 ; k2 < _n ; ++k2 )
157 {
158 if ( !_fixed_vars[k1] && !_fixed_vars[k2] )
159 _index[c2++] = c1;
160 ++c1;
161 }
162 }
163
164 /*---------------------------------------------------------*/
165 /* check evaluation point outputs before the integration */
166 /* into an interpolation set (private) */
167 /*---------------------------------------------------------*/
check_outputs(const NOMAD::Point & bbo,int m) const168 bool NOMAD::Quad_Model::check_outputs ( const NOMAD::Point & bbo , int m ) const {
169
170 if ( bbo.size() != m )
171 return false;
172
173 for ( int i = 0 ; i < m ; ++i )
174 if ( !bbo[i].is_defined() || bbo[i].value() > NOMAD::MODEL_MAX_OUTPUT )
175 return false;
176
177 return true;
178 }
179
180 /*-----------------------------------------------------------*/
181 /* construct the interpolation set Y */
182 /*-----------------------------------------------------------*/
construct_Y(const NOMAD::Point & center,const NOMAD::Point & interpolation_radius,int max_Y_size)183 void NOMAD::Quad_Model::construct_Y ( const NOMAD::Point & center ,
184 const NOMAD::Point & interpolation_radius ,
185 int max_Y_size )
186 {
187 _error_flag = true;
188
189 if ( center.size() != _n ||
190 interpolation_radius.size() != _n ||
191 !center.is_complete() ||
192 !interpolation_radius.is_complete() )
193 return;
194
195 _error_flag = false;
196 _center = center;
197
198 int m = static_cast<int> ( _bbot.size() );
199
200 // browse the cache:
201 const NOMAD::Eval_Point * cur = _cache.begin();
202 while ( cur )
203 {
204
205 if ( cur->get_eval_status() == NOMAD::EVAL_OK &&
206 cur->get_n () == _n &&
207 _signature.is_compatible (*cur) )
208 {
209
210 const NOMAD::Point & bbo = cur->get_bb_outputs();
211
212 if ( check_outputs ( bbo , m ) )
213 {
214
215 // the center point has been found
216 // (it is put in first position):
217 if ( _center == *cur )
218 {
219 _Y.push_back ( new NOMAD::Eval_Point ( *cur ) );
220 int nYm1 = get_nY() - 1;
221 if ( nYm1 > 0 )
222 {
223 NOMAD::Eval_Point * tmp = _Y[0];
224 _Y[0 ] = _Y[nYm1];
225 _Y[nYm1] = tmp;
226 }
227 }
228
229 // other points must within the interpolation radius:
230 else if ( is_within_radius ( *cur , interpolation_radius ) )
231 {
232 _Y.push_back ( new NOMAD::Eval_Point ( *cur ) );
233 }
234 }
235 }
236 cur = _cache.next();
237 }
238
239 // respect the limit on the number of points:
240 if ( get_nY() > max_Y_size )
241 reduce_Y ( center , max_Y_size );
242 }
243
244 /*-----------------------------------------------------------------*/
245 /* reduce the number of interpolation points */
246 /*-----------------------------------------------------------------*/
reduce_Y(const NOMAD::Point & center,int max_Y_size)247 void NOMAD::Quad_Model::reduce_Y ( const NOMAD::Point & center ,
248 int max_Y_size )
249 {
250 int nY = get_nY();
251
252 if ( nY <= max_Y_size )
253 return;
254
255 std::multiset<NOMAD::Model_Sorted_Point> Ys;
256 for ( int k = 0 ; k < nY ; ++k )
257 Ys.insert ( NOMAD::Model_Sorted_Point ( _Y[k] , center ) );
258
259 _Y.clear();
260
261 std::multiset<NOMAD::Model_Sorted_Point>::const_iterator it , end = Ys.end();
262 for ( it = Ys.begin() ; it != end ; ++it ) {
263 if ( get_nY() < max_Y_size )
264 _Y.push_back ( static_cast<NOMAD::Eval_Point *> ( it->get_point() ) );
265 else
266 delete it->get_point();
267 }
268 }
269
270 /*-----------------------------------------------------------*/
271 /* check if an unscaled point is in B(center,radius) for a */
272 /* given radius (private) */
273 /*-----------------------------------------------------------*/
is_within_radius(const NOMAD::Point & x,const NOMAD::Point & radius) const274 bool NOMAD::Quad_Model::is_within_radius ( const NOMAD::Point & x ,
275 const NOMAD::Point & radius ) const
276 {
277 if ( x.size() != _n || radius.size() != _n )
278 return false;
279
280 for ( int i = 0 ; i < _n ; ++i )
281 {
282 if ( !x[i].is_defined() ||
283 !radius[i].is_defined() ||
284 radius[i] < ( x[i] - _center[i]).abs() )
285 return false;
286 }
287 return true;
288 }
289
290 /*------------------------------------------------------*/
291 /* check if a scaled point is inside the trust radius */
292 /*------------------------------------------------------*/
is_within_trust_radius(const NOMAD::Point & x) const293 bool NOMAD::Quad_Model::is_within_trust_radius ( const NOMAD::Point & x ) const
294 {
295 // check that all scaled coordinates are in [-1;1] and
296 // that fixed variables are equal to zero:
297 for ( int i = 0 ; i < _n ; ++i )
298 if ( !_ref [i].is_defined() ||
299 !_scaling[i].is_defined() ||
300 ! x[i].is_defined() ||
301 x[i].abs() > 1.0 ||
302 ( _fixed_vars[i] && x[i] != 0.0 ) )
303 return false;
304 return true;
305 }
306
307 /*--------------------------------------------------------------*/
308 /* . define scaling to put all coordinates centered in [-r;r] */
309 /* . looks also for fixed variables */
310 /*--------------------------------------------------------------*/
define_scaling(const NOMAD::Double & r)311 void NOMAD::Quad_Model::define_scaling ( const NOMAD::Double & r )
312 {
313 if ( _error_flag || _Y.empty() ) {
314 _error_flag = true;
315 return;
316 }
317
318 int i , k;
319 int nY = get_nY();
320 NOMAD::Point min(_n) , max(_n);
321 NOMAD::Double tmp;
322
323 // The parameters defining the scaling with rotation (see define_scaling_by_direction) are cleared.
324 // Only the parameters for the basic scaling are set
325 _dirP.clear();
326 _epsilon.clear();
327 _delta_m.clear();
328
329
330 _scaling.clear();
331 _ref.clear ();
332 _ref.reset ( _n );
333 _scaling.reset ( _n );
334
335 // compute the reference (center of Y):
336 for ( k = 0 ; k < nY ; ++k )
337 {
338
339 if ( !_Y[k] || _n != _Y[k]->size() )
340 {
341 _error_flag = true;
342 return;
343 }
344
345 for ( i = 0 ; i < _n ; ++i )
346 {
347 tmp = (*_Y[k])[i];
348 if ( !min[i].is_defined() || tmp < min[i] )
349 min[i] = tmp;
350 if ( !max[i].is_defined() || tmp > max[i] )
351 max[i] = tmp;
352 }
353 }
354
355 for ( i = 0 ; i < _n ; ++i )
356 _ref[i] = ( min[i] + max[i] ) / 2.0;
357
358 #ifdef MODEL_STATS
359 _Yw = NOMAD::Double();
360 for ( i = 0 ; i < _n ; ++i )
361 {
362 tmp = max[i]-min[i];
363 if ( !_Yw.is_defined() || tmp > _Yw )
364 _Yw = tmp;
365 }
366 #endif
367
368 #ifdef DEBUG
369 _out << std::endl
370 << "define_scaling(): reference = ( " << _ref << " )" << std::endl;
371 #endif
372
373 // compute the scaling (and detect fixed variables):
374 for ( k = 0 ; k < nY ; ++k ) {
375
376 for ( i = 0 ; i < _n ; ++i ) {
377 tmp = ( (*_Y[k])[i] - _ref[i] ).abs();
378 if ( !_scaling[i].is_defined() || _scaling[i] < tmp )
379 _scaling[i] = tmp;
380 }
381 }
382
383 _nfree = _n;
384
385 for ( i = 0 ; i < _n ; ++i )
386 {
387 if ( _scaling[i] == 0.0 )
388 {
389 _scaling [i] = 0.0;
390 _fixed_vars[i] = true;
391 --_nfree;
392 if ( _nfree == 0 )
393 {
394 _scaling.clear();
395 _ref.clear();
396 _error_flag = true;
397 return;
398 }
399 }
400 else
401 _scaling[i] *= 1.0/r; // all coordinates in [-r;r]
402 }
403
404 if ( _nfree < _n )
405 init_alpha();
406
407 for ( k = 0 ; k < nY ; ++k )
408 {
409 if ( !scale ( *_Y[k] ) )
410 {
411 _scaling.clear();
412 _error_flag = true;
413 return;
414 }
415 }
416
417 #ifdef DEBUG
418 _out << "define_scaling(): scaling = ( " << _scaling << " )" << std::endl;
419 #endif
420
421 _error_flag = false;
422 }
423
424 /*-------------------------------------------------------------------*/
425 /* . Scaling with rotation based on a set of directions. */
426 /* See paper: */
427 /* Reducing the number of function evaluations in */
428 /* Mesh Adaptive Direct Search algorithms, Audet, Ianni, */
429 /* LeDigabel, Tribes, 2014 */
430 /* . looks also for fixed variables */
431 /*-------------------------------------------------------------------*/
define_scaling_by_directions(const std::list<NOMAD::Direction> & dirs,const NOMAD::Point & delta_m,const NOMAD::Double & epsilon)432 void NOMAD::Quad_Model::define_scaling_by_directions ( const std::list<NOMAD::Direction> & dirs, const NOMAD::Point & delta_m, const NOMAD::Double & epsilon )
433 {
434 if ( _error_flag || _Y.empty() ) {
435 _error_flag = true;
436 return;
437 }
438
439 int i , k;
440 int nY = get_nY();
441 NOMAD::Point min(_n) , max(_n);
442 NOMAD::Double tmp;
443
444
445 // The parameters defining the basic scaling (see define_scaling) are cleared.
446 // Only the parameters for the direction scaling are set
447 _scaling.clear();
448 _ref.clear ();
449
450 // For direction scaling
451 if (static_cast<int> (dirs.size())!=_n || static_cast<int>(delta_m.size()) != _n || epsilon<=0.0 || epsilon>=1)
452 {
453 _error_flag = true;
454 return;
455 }
456 _delta_m=delta_m;
457 // Get D' from dirs (scaling with delta_m
458 std::list<NOMAD::Direction>::const_iterator itDir;
459 for (itDir=dirs.begin(); itDir != dirs.end(); itDir++)
460 {
461 NOMAD::Direction dir_i(_n,0.0,itDir->get_type());
462 dir_i.set_index(itDir->get_index());
463 for ( int i = 0 ; i < _n ; ++i )
464 {
465 if (_delta_m[i]<=0.0)
466 {
467 _error_flag=true;
468 return;
469 }
470 dir_i[i]=(*itDir)[i]/_delta_m[i];
471 }
472 _dirP.push_back(dir_i);
473 }
474
475 _epsilon=epsilon;
476
477
478 // compute the min and the max:
479 for ( k = 0 ; k < nY ; ++k )
480 {
481
482 if ( !_Y[k] || _n != _Y[k]->size() )
483 {
484 _error_flag = true;
485 return;
486 }
487
488 for ( i = 0 ; i < _n ; ++i )
489 {
490 tmp = (*_Y[k])[i];
491 if ( !min[i].is_defined() || tmp < min[i] )
492 min[i] = tmp;
493 if ( !max[i].is_defined() || tmp > max[i] )
494 max[i] = tmp;
495 }
496 }
497
498 #ifdef MODEL_STATS
499 _Yw = NOMAD::Double();
500 for ( i = 0 ; i < _n ; ++i )
501 {
502 tmp = max[i]-min[i];
503 if ( !_Yw.is_defined() || tmp > _Yw )
504 _Yw = tmp;
505 }
506 #endif
507
508 // Detect fixed variables:
509 _nfree = _n;
510 for ( i = 0 ; i < _n ; ++i )
511 {
512 bool fixed_var_i=true;
513 for ( k = 0 ; k < nY ; ++k )
514 {
515 if ( ( (*_Y[k])[i] - _center[i] ).abs() > 0.0 )
516 {
517 fixed_var_i=false;
518 break;
519 }
520 }
521 _fixed_vars[i]=fixed_var_i;
522 if (fixed_var_i)
523 --_nfree;
524
525 if ( _nfree == 0 )
526 {
527 _scaling.clear();
528 _ref.clear();
529 _dirP.clear();
530 _error_flag = true;
531 return;
532 }
533 }
534 if ( _nfree < _n )
535 init_alpha();
536
537 // Perform scaling of Y
538 for ( k = 0 ; k < nY ; ++k )
539 {
540 if ( !scale ( *_Y[k] ) )
541 {
542 _scaling.clear();
543 _dirP.clear();
544 _error_flag = true;
545 return;
546 }
547 }
548
549 #ifdef DEBUG
550 _out << "define_scaling_by_direction(): " << std::endl;
551 for ( itDir = _dirP.begin() ; itDir != _dirP.end() ; ++itDir )
552 {
553 _out << "dirPrime ";
554 _out.display_int_w ( (*itDir).get_index() , static_cast<int>(_dirP.size()) );
555 _out << " : " << *itDir << std::endl;
556 }
557 #endif
558
559 _error_flag = false;
560 }
561
562
563 /*--------------------------------------------------------------*/
564 /* scale a point */
565 /*--------------------------------------------------------------*/
scale(NOMAD::Point & x) const566 bool NOMAD::Quad_Model::scale ( NOMAD::Point & x ) const
567 {
568 if ( _error_flag || _n != x.size() )
569 return false;
570
571 if (_dirP.size()==0)
572 {
573 // Scale without rotation
574 for ( int i = 0 ; i < _n ; ++i )
575 {
576 if ( !_ref [i].is_defined() ||
577 !_scaling[i].is_defined() ||
578 ! x[i].is_defined() )
579 return false;
580 x[i] -= _ref[i];
581 if ( _scaling[i] != 0 )
582 x[i] /= _scaling[i];
583 }
584 }
585 else
586 {
587 if (! _epsilon.is_defined() || !_delta_m.is_complete())
588 return false;
589 // Scale with rotation based on direction and center (see paper Reducing the number of function evaluations in Mesh Adaptive Direct Search algorithms, Audet, Ianni, LeDigabel, Tribes, 2014
590 // T(y)=(D')^-1*(center-x)/delta_m/(1-epsilon) - epsilon*1/(1-epsilon)
591 // (D')^-1=(D')^T/normCol^2
592 NOMAD::Point temp(_n,0.0);
593 NOMAD::Double normCol2=0.0;
594 std::list<NOMAD::Direction>::const_iterator itDir=_dirP.begin();
595 for ( int i = 0 ; i < _n ; ++i )
596 {
597 normCol2+=pow((*itDir)[i].value(),2.0);
598
599 if (_delta_m[i] !=0.0)
600 temp[i]=(_center[i].value()-x[i].value())/_delta_m[i].value()/(1.0-_epsilon.value());
601 else
602 return false;
603 x[i]=0.0;
604 }
605 int j=0;
606 for (itDir=_dirP.begin(); itDir != _dirP.end(); itDir++,j++)
607 {
608 for ( int i = 0 ; i < _n ; ++i )
609 {
610 x[j]+=temp[i].value()*(*itDir)[i].value()/normCol2.value();
611 }
612 x[j]-=_epsilon.value()/(1.0-_epsilon.value());
613 }
614 }
615
616 return true;
617 }
618
619 /*-----------------------------------------------------------*/
620 /* unscale a point */
621 /*-----------------------------------------------------------*/
unscale(NOMAD::Point & x) const622 bool NOMAD::Quad_Model::unscale ( NOMAD::Point & x ) const
623 {
624 if ( _error_flag || _n != x.size() )
625 return false;
626
627 if (_dirP.size()==0)
628 {
629 // Scale without rotation
630 for ( int i = 0 ; i < _n ; ++i )
631 {
632 if ( !_ref [i].is_defined() ||
633 !_scaling[i].is_defined() ||
634 ! x[i].is_defined() )
635 return false;
636
637 x[i] *= _scaling[i];
638 x[i] += _ref [i];
639 }
640 }
641 else
642 {
643
644 if (! _epsilon.is_defined() || !_delta_m.is_complete())
645 return false;
646
647 // UnScale with rotation see paper Reducing the number of function evaluations in Mesh Adaptive Direct Search algorithms, Audet, Ianni, LeDigabel, Tribes, 2014
648 //T^−1(x) = center+ _delta_m Dp ((ε−1)x−ε1)
649 NOMAD::Point temp(_n,0.0);
650 for ( int i = 0 ; i < _n ; ++i )
651 {
652 temp[i]=(x[i]*(_epsilon-1.0)-_epsilon)*_delta_m[i];
653 x[i]=0.0;
654 }
655 std::list<NOMAD::Direction>::const_iterator itDir;
656 int j=0;
657 for (itDir=_dirP.begin(); itDir != _dirP.end(); itDir++,j++)
658 {
659 for (int i=0 ; i< _n ; i++)
660 {
661 x[i]+=temp[j]*(*itDir)[i];
662 }
663 }
664 for ( int i = 0 ; i < _n ; ++i )
665 {
666 x[i]+=_center[i];
667 }
668 }
669
670 return true;
671 }
672
673 /*-----------------------------------------------------------*/
674 /* unscale the slope at a point */
675 /*-----------------------------------------------------------*/
unscale_grad(NOMAD::Point & x) const676 bool NOMAD::Quad_Model::unscale_grad ( NOMAD::Point & x ) const
677 {
678 if ( _error_flag || _n != x.size() )
679 return false;
680
681 for ( int i = 0 ; i < _n ; ++i )
682 {
683
684 if (!_scaling[i].is_defined() || !x[i].is_defined() )
685 return false;
686
687 x[i] *= _scaling[i];
688 }
689
690 return true;
691 }
692
693 /*------------------------------------------------------------------*/
694 /* compute the element (i,j) of the interpolation matrix M(phi,Y) */
695 /* (private) */
696 /*------------------------------------------------------------------*/
compute_M(int i,int j) const697 double NOMAD::Quad_Model::compute_M ( int i , int j ) const {
698
699 if ( _error_flag )
700 return 0.0;
701
702 if ( j == 0 )
703 return 1.0;
704
705 if ( j <= _nfree )
706 return (*_Y[i])[_index[j]-1].value();
707
708 if ( j <= 2 * _nfree )
709 return 0.5 * pow ( (*_Y[i])[_index[j-_nfree]-1].value() , 2.0 );
710
711 int nm1 = _nfree - 1;
712 int jm2n , dec , r , i1 , i2;
713
714 jm2n = j - 2 * _nfree;
715 dec = nm1;
716 r = jm2n;
717 i1 = -1;
718
719 while ( r > 0 ) {
720 r -= dec;
721 ++i1;
722 --dec;
723 }
724
725 i2 = r + nm1;
726
727 return (*_Y[i])[_index[i1+1]-1].value() * (*_Y[i])[_index[i2+1]-1].value();
728 }
729
730 /*-----------------------------------------------------------*/
731 /* construct m models (one by output) */
732 /*-----------------------------------------------------------*/
construct(bool use_WP,double eps,int max_mpn,int max_Y_size)733 void NOMAD::Quad_Model::construct ( bool use_WP ,
734 double eps ,
735 int max_mpn ,
736 int max_Y_size )
737 {
738 if ( _error_flag )
739 return;
740
741 int p1 = get_nY();
742
743
744 // MFN interpolation:
745 if ( p1 < _n_alpha ) {
746 _interpolation_type = NOMAD::MFN;
747 _error_flag = !construct_MFN_model ( eps , max_mpn , max_Y_size );
748 }
749 else {
750
751 _error_flag = true;
752
753 // well-poised regression:
754 if ( use_WP && p1 > _n_alpha ) {
755 _interpolation_type = NOMAD::WP_REGRESSION;
756 _error_flag = !construct_WP_model ( max_Y_size );
757 }
758
759 // regression:
760 if ( _error_flag ) {
761 _interpolation_type = NOMAD::REGRESSION;
762 _error_flag = !construct_regression_model ( eps , max_mpn , max_Y_size );
763 }
764 }
765 }
766
767 /*---------------------------------------------------------------*/
768 /* find interpolation point with max Lagrange polynomial value */
769 /*---------------------------------------------------------------*/
770 /* . ji = argmax |li(x)| for x in Y */
771 /* . used in construct_WP_model() */
772 /* . private */
773 /*---------------------------------------------------------------*/
find_max_lix(const NOMAD::Point & li,const std::vector<NOMAD::Eval_Point * > & Y,int i1,int i2,NOMAD::Double & max_lix) const774 int NOMAD::Quad_Model::find_max_lix
775 ( const NOMAD::Point & li ,
776 const std::vector<NOMAD::Eval_Point *> & Y ,
777 int i1 ,
778 int i2 ,
779 NOMAD::Double & max_lix ) const
780 {
781 max_lix = -1.0;
782 int ji = -1;
783 NOMAD::Double tmp;
784 for ( int j = i1 ; j <= i2 ; ++j ) {
785 tmp = eval ( *Y[j] , li );
786 if ( tmp.is_defined() ) {
787 tmp = tmp.abs();
788 if ( tmp > max_lix ) {
789 max_lix = tmp;
790 ji = j;
791 }
792 }
793 }
794 if ( ji < 0 )
795 max_lix.clear();
796 return ji;
797 }
798
799 /*-----------------------------------------------------------*/
800 /* construct well-poised (WP) model (private) */
801 /*-----------------------------------------------------------*/
construct_WP_model(int max_Y_size)802 bool NOMAD::Quad_Model::construct_WP_model ( int max_Y_size )
803 {
804
805 #ifdef DEBUG
806 _out << std::endl
807 << NOMAD::open_block ( "NOMAD::Quad_Model::construct_WP_model()" );
808 #endif
809
810 // check the set Y:
811 if ( !check_Y() )
812 return false;
813
814 int i , j , k , p1 = get_nY();
815
816 // the number of points (p+1) must be in [1+(n+1)(n+2)/2;MS_MAX_Y_SIZE]:
817 if ( p1 <= _n_alpha || p1 > max_Y_size ) {
818 #ifdef DEBUG
819 _out << std::endl
820 << "NOMAD::Quad_Model::construct_WP_model(): "
821 << "(p+1) not in [1+(n+1)(n+2)/2;" << max_Y_size << "]"
822 << std::endl << NOMAD::close_block() << std::endl;
823 #endif
824 return false;
825 }
826
827 // Lagrange polynomials:
828 std::vector<NOMAD::Point *> l;
829 for ( i = 0 ; i < _n_alpha ; ++i ) {
830 l.push_back ( new NOMAD::Point ( _n_alpha ) );
831 for ( j = 0 ; j < _n_alpha ; ++j )
832 (*l[i])[j] = (i==j) ? 1.0 : 0.0;
833 }
834
835 // creation of sets Y1 and Y2; Y2 contains all available points
836 // of _Y and Y1 will be the 'well-poised' set with n_alpha points:
837 std::vector<NOMAD::Eval_Point *> Y1 , Y2 = _Y;
838 int iy2 , ny2m1 = p1-1;
839 NOMAD::Double max_lix , liyi , ljyi;
840
841 // we init Y1 with the first point of Y:
842 Y1.push_back ( Y2[0] );
843 Y2[0] = Y2[ny2m1];
844 Y2.resize ( ny2m1 );
845
846 // use algo 6.2 p.95 of the DFO book in order to construct Lagrange polynomials:
847 // -----------------------------------------------------------------------------
848 for ( i = 0 ; i < _n_alpha ; ++i ) {
849
850 // 1. point selection (select a point in Y2: Y2[iy2]):
851 // -------------------
852 if ( i > 0 ) {
853
854 ny2m1 = static_cast<int>(Y2.size())-1;
855 iy2 = find_max_lix ( *l[i] , Y2 , 0 , ny2m1 , max_lix );
856
857 if ( iy2 < 0 ) {
858 #ifdef DEBUG
859 _out << std::endl
860 << "NOMAD::Quad_Model::construct_WP_model(): "
861 << "cannot find candidate in Y"
862 << std::endl << NOMAD::close_block() << std::endl;
863 #endif
864 for ( i = 0 ; i < _n_alpha ; ++i )
865 delete l[i];
866 return false;
867 }
868
869 // add Y2[iy2] in Y1:
870 Y1.push_back ( Y2[iy2] );
871 Y2[iy2] = Y2[ny2m1];
872 Y2.resize (ny2m1);
873 }
874
875 // 2. normalization:
876 // -----------------
877 liyi = eval ( *Y1[i] , *l[i] );
878
879 if ( liyi.abs().value() < 1e-15 ) {
880 #ifdef DEBUG
881 _out << std::endl
882 << "NOMAD::Quad_Model::construct_WP_model(): set Y is not poised"
883 << std::endl << NOMAD::close_block() << std::endl;
884 #endif
885 for ( i = 0 ; i < _n_alpha ; ++i )
886 delete l[i];
887 return false;
888 }
889
890 for ( k = 0 ; k < _n_alpha ; ++k ) {
891 (*l[i])[k] /= liyi;
892 if ( (*l[i])[k].abs().value() < 1e-15 )
893 (*l[i])[k] = 0.0;
894 }
895
896 // 3. orthogonalization:
897 // ---------------------
898 for ( j = 0 ; j < _n_alpha ; ++j )
899 if ( j != i ) {
900 ljyi = eval ( *Y1[i] , *l[j] );
901 for ( k = 0 ; k < _n_alpha ; ++k ) {
902 (*l[j])[k] = (*l[j])[k] - ljyi * (*l[i])[k];
903 if ( (*l[j])[k].abs().value() < 1e-15 )
904 (*l[j])[k] = 0.0;
905 }
906 }
907 }
908
909 #ifdef DEBUG
910 display_lagrange_polynomials ( l , Y1 );
911 #endif
912
913 // compute alpha:
914 // --------------
915 int m = static_cast<int> ( _bbot.size() );
916 for ( i = 0 ; i < m ; ++i )
917 if ( _alpha[i] ) {
918 for ( j = 0 ; j < _n_alpha ; ++j ) {
919 (*_alpha[i])[j] = 0.0;
920 for ( k = 0 ; k < _n_alpha ; ++k )
921 (*_alpha[i])[j] += Y1[k]->get_bb_outputs()[i] * (*l[k])[j];
922 }
923 }
924
925 // poisedness improvement using algorithm 6.3 page 95:
926 // ---------------------------------------------------
927
928 // old alpha:
929 NOMAD::Point ** old_alpha = new NOMAD::Point * [m] , ** tmp_alpha;
930 for ( i = 0 ; i < m ; ++i )
931 old_alpha[i] = ( _alpha[i] ) ?
932 new NOMAD::Point ( _n_alpha ) : NULL;
933
934 int ik;
935 NOMAD::Double ljyk , lkyk , lix , new_rel_err ,
936 cur_rel_err = compute_max_rel_err();
937
938 if ( cur_rel_err.is_defined() && cur_rel_err.value() > 1e-15 ) {
939
940 for ( int niter = 0 ; niter < 10 ; ++niter ) {
941
942 ny2m1 = static_cast<int>(Y2.size())-1;
943
944 if ( ny2m1 < 0 )
945 break;
946
947 max_lix = -1.0;
948 iy2 = -1;
949 ik = -1;
950
951 for ( i = 0 ; i < _n_alpha ; ++i ) {
952
953 j = find_max_lix ( *l[i] , Y2 , 0 , ny2m1 , lix );
954 if ( j >= 0 && lix > max_lix ) {
955 max_lix = lix;
956 iy2 = j;
957 ik = i;
958 }
959 }
960
961 if ( ik < 0 )
962 break;
963
964 // set Y1[ik] = Y2[iy2]:
965 Y1[ik ] = Y2[iy2];
966 Y2[iy2] = Y2[ny2m1];
967 Y2.resize ( ny2m1 );
968
969 lkyk = eval ( *Y1[ik] , *l[ik] );
970
971 if ( lkyk.abs() <= 1e-15 )
972 break;
973
974 // update Lagrange polynomials:
975 // ----------------------------
976
977 // normalization and orthogonalization:
978 for ( i = 0 ; i < _n_alpha ; ++i )
979 (*l[ik])[i] /= lkyk;
980
981 for ( j = 0 ; j < _n_alpha ; ++j ) {
982 if ( j != ik ) {
983 ljyk = eval ( *Y1[ik] , *l[j] );
984 for ( i = 0 ; i < _n_alpha ; ++i )
985 (*l[j])[i] = (*l[j])[i] - ljyk * (*l[ik])[i];
986 }
987 }
988
989 // save old alpha and compute new one:
990 for ( i = 0 ; i < m ; ++i )
991 if ( _alpha[i] ) {
992 *(old_alpha[i]) = *(_alpha[i]);
993 for ( j = 0 ; j < _n_alpha ; ++j ) {
994 (*_alpha[i])[j] = 0.0;
995 for ( k = 0 ; k < _n_alpha ; ++k )
996 (*_alpha[i])[j] += Y1[k]->get_bb_outputs()[i] * (*l[k])[j];
997 }
998 }
999
1000 // compute new error:
1001 new_rel_err = compute_max_rel_err();
1002
1003 // if no better error, restore old alpha and exit loop:
1004 if ( !new_rel_err.is_defined() || new_rel_err >= cur_rel_err ) {
1005 tmp_alpha = _alpha;
1006 _alpha = old_alpha;
1007 old_alpha = tmp_alpha;
1008 break;
1009 }
1010
1011 cur_rel_err = new_rel_err;
1012 }
1013 }
1014
1015 for ( i = 0 ; i < m ; ++i )
1016 delete old_alpha[i];
1017 delete [] old_alpha;
1018 for ( i = 0 ; i < _n_alpha ; ++i )
1019 delete l[i];
1020
1021 #ifdef DEBUG
1022 _out.close_block();
1023 #endif
1024
1025 return true;
1026 }
1027
1028 /*-----------------------------------------------------------*/
1029 /* construct regression model (private) */
1030 /*-----------------------------------------------------------*/
construct_regression_model(double eps,int max_mpn,int max_Y_size)1031 bool NOMAD::Quad_Model::construct_regression_model ( double eps ,
1032 int max_mpn ,
1033 int max_Y_size )
1034 {
1035 #ifdef DEBUG
1036 _out << std::endl
1037 << NOMAD::open_block ( "NOMAD::Quad_Model::construct_regression_model()" );
1038 #endif
1039
1040 _error_flag = false;
1041
1042 // check the set Y:
1043 if ( !check_Y() )
1044 return false;
1045
1046 int p1 = get_nY();
1047
1048 // the number of points (p+1) must be in [(n+1)(n+2)/2;MS_MAX_Y_SIZE]:
1049 if ( p1 < _n_alpha || p1 > max_Y_size ) {
1050 #ifdef DEBUG
1051 _out << std::endl
1052 << "NOMAD::Quad_Model::construct_regression_model(): "
1053 << "(p+1) not in [(n+1)(n+2)/2;"
1054 << max_Y_size << "]"
1055 << std::endl << NOMAD::close_block() << std::endl;
1056 #endif
1057 return false;
1058 }
1059
1060 // for this procedure, the number of points is limited to 500
1061 // (because of the SVD decomposition):
1062 if ( p1 > 500 ) {
1063 reduce_Y ( NOMAD::Point ( _n , 0.0 ) , 500 );
1064 p1 = 500;
1065 }
1066
1067 // construct the matrix F=M'M (_n_alpha,_n_alpha):
1068 // -----------------------------------------------
1069 int i , j , k;
1070 double ** F = new double *[_n_alpha];
1071 double ** M = new double *[p1];
1072 for ( i = 0 ; i < p1 ; ++i ) {
1073 M[i] = new double[_n_alpha];
1074 for ( j = 0 ; j < _n_alpha ; ++j )
1075 M[i][j] = compute_M ( i , j );
1076 }
1077
1078 for ( i = 0 ; i < _n_alpha ; ++i ) {
1079 F[i] = new double[_n_alpha];
1080 for ( j = 0 ; j <= i ; ++j ) {
1081 F[i][j] = 0.0;
1082 for ( k = 0 ; k < p1 ; ++k )
1083 F[i][j] += M[k][i] * M[k][j];
1084 if ( i != j )
1085 F[j][i] = F[i][j];
1086 }
1087 }
1088
1089 #ifdef DEBUG
1090 _out << std::endl << "F=";
1091 for ( i = 0 ; i < _n_alpha ; ++i ) {
1092 _out << "\t";
1093 for ( j = 0 ; j < _n_alpha ; ++j )
1094 _out << std::setw(12) << F[i][j] << " ";
1095 _out << std::endl;
1096 }
1097 #endif
1098
1099 bool error = false;
1100
1101 // SVD decomposition of the F matrix (F=U.W.V'):
1102 // ---------------------------------------------
1103 // (F will be transformed in U)
1104
1105 double * W = new double [_n_alpha];
1106 double ** V = new double *[_n_alpha];
1107 for ( i = 0 ; i < _n_alpha ; ++i )
1108 V[i] = new double[_n_alpha];
1109
1110 std::string error_msg;
1111 if ( NOMAD::SVD_decomposition ( error_msg , F , W , V , _n_alpha , _n_alpha , max_mpn ) ) {
1112
1113 // compute condition number:
1114 compute_cond ( W , _n_alpha , eps );
1115
1116 #ifdef DEBUG
1117 _out << std::endl << "F=";
1118 for ( i = 0 ; i < _n_alpha ; ++i ) {
1119 _out << "\t";
1120 for ( j = 0 ; j < _n_alpha ; ++j )
1121 _out << std::setw(12) << F[i][j] << " ";
1122 _out << std::endl;
1123 }
1124
1125 _out << std::endl << "W=\t";
1126 for ( i = 0 ; i < _n_alpha ; ++i )
1127 _out << std::setw(12) << W[i] << " ";
1128 _out << std::endl << std::endl << "cond=" << _cond << std::endl;
1129
1130 _out << std::endl << "V=";
1131 for ( i = 0 ; i < _n_alpha ; ++i ) {
1132 _out << "\t";
1133 for ( j = 0 ; j < _n_alpha ; ++j )
1134 _out << std::setw(12) << V[i][j] << " ";
1135 _out << std::endl;
1136 }
1137 #endif
1138
1139 }
1140 else {
1141 #ifdef DEBUG
1142 _out << std::endl << "NOMAD::Quad_Model::construct_regression_model(): "
1143 << "SVD decomposition (" << error_msg << ")"
1144 << std::endl << NOMAD::close_block() << std::endl;
1145 #endif
1146 error = true;
1147 _cond.clear();
1148 }
1149
1150 // resolution of system F.alpha = M'.f(Y):
1151 // ---------------------------------------
1152 if ( !error ) {
1153 int m = static_cast<int> ( _bbot.size() );
1154 for ( i = 0 ; i < m ; ++i )
1155 if ( _alpha[i] )
1156 solve_regression_system ( M , F , W , V , i , *_alpha[i] , eps );
1157 }
1158
1159 // free memory:
1160 for ( i = 0 ; i < _n_alpha ; ++i ) {
1161 delete [] F[i];
1162 delete [] V[i];
1163 }
1164 for ( i = 0 ; i < p1 ; ++i )
1165 delete [] M[i];
1166 delete [] M;
1167 delete [] F;
1168 delete [] V;
1169 delete [] W;
1170
1171 #ifdef DEBUG
1172 _out.close_block();
1173 #endif
1174
1175 return !error;
1176 }
1177
1178 /*-------------------------------------------------------------*/
1179 /* compute condition number (private) */
1180 /*-------------------------------------------------------------*/
compute_cond(const double * W,int n,double eps)1181 void NOMAD::Quad_Model::compute_cond ( const double * W , int n , double eps )
1182 {
1183 double min = NOMAD::INF;
1184 double max = -min;
1185 for ( int i = 0 ; i < n ; ++i ) {
1186 if ( W[i] < min )
1187 min = W[i];
1188 if ( W[i] > max )
1189 max = W[i];
1190 }
1191 if ( min < eps )
1192 min = eps;
1193 _cond = max / min;
1194 }
1195
1196 /*-------------------------------------------------------------*/
1197 /* resolution of system F.alpha = M'.f(Y) for the regression */
1198 /* (private) */
1199 /*-------------------------------------------------------------*/
solve_regression_system(double ** M,double ** F,double * W,double ** V,int bbo_index,NOMAD::Point & alpha,double eps) const1200 void NOMAD::Quad_Model::solve_regression_system ( double ** M ,
1201 double ** F ,
1202 double * W ,
1203 double ** V ,
1204 int bbo_index ,
1205 NOMAD::Point & alpha ,
1206 double eps ) const
1207 {
1208 // resize the alpha vector:
1209 if ( alpha.size() != _n_alpha )
1210 alpha.reset ( _n_alpha , 0.0 );
1211
1212 double * alpha_tmp = new double [_n_alpha];
1213 int i , k , p1 = get_nY();
1214
1215 // solve the system:
1216 for ( i = 0 ; i < _n_alpha ; ++i ) {
1217 alpha_tmp[i] = 0.0;
1218 for ( k = 0 ; k < p1 ; ++k )
1219 alpha_tmp[i] += M[k][i] * ( _Y[k]->get_bb_outputs()[bbo_index].value() );
1220 }
1221
1222 double * alpha_tmp2 = new double [_n_alpha];
1223
1224 // some W values will be zero (or near zero);
1225 // each value that is smaller than eps is ignored
1226
1227 for ( i = 0 ; i < _n_alpha ; ++i ) {
1228 alpha_tmp2[i] = 0.0;
1229 for ( k = 0 ; k < _n_alpha ; ++k )
1230 if ( W[i] > eps )
1231 alpha_tmp2[i] += F[k][i] * alpha_tmp[k] / W[i];
1232 }
1233
1234 delete [] alpha_tmp;
1235
1236 for ( i = 0 ; i < _n_alpha ; ++i ) {
1237 alpha[i] = 0.0;
1238 for ( k = 0 ; k < _n_alpha ; ++k )
1239 alpha[i] += V[i][k] * alpha_tmp2[k];
1240 }
1241
1242 delete [] alpha_tmp2;
1243 }
1244
1245 /*----------------------------------------------------------*/
1246 /* construct Minimum Frobenius Norm (MFN) model (private) */
1247 /*----------------------------------------------------------*/
construct_MFN_model(double eps,int max_mpn,int max_Y_size)1248 bool NOMAD::Quad_Model::construct_MFN_model ( double eps ,
1249 int max_mpn ,
1250 int max_Y_size )
1251 {
1252 #ifdef DEBUG
1253 _out << std::endl
1254 << NOMAD::open_block ( "NOMAD::Quad_Model::construct_MFN_model()" );
1255 #endif
1256
1257 // check the set Y:
1258 if ( !check_Y() )
1259 return false;
1260
1261 int p1 = get_nY();
1262
1263 // the number of points (p+1) must be in [n+1;(n+1)(n+2)/2-1]:
1264 if ( p1 <= _nfree || p1 >= _n_alpha ) {
1265 #ifdef DEBUG
1266 _out << std::endl
1267 << "NOMAD::Quad_Model::construct_MFN_model(): "
1268 << "(p+1) not in [n+1;(n+1)(n+2)/2-1]"
1269 << std::endl << NOMAD::close_block() << std::endl;
1270 #endif
1271 return false;
1272 }
1273
1274 // for this procedure, the number of points is limited to 250
1275 // (because of the SVD decomposition):
1276 if ( p1 > 250 ) {
1277 reduce_Y ( NOMAD::Point ( _n , 0.0 ) , 250 );
1278 p1 = 250;
1279 }
1280
1281 // construct the matrix F (4 parts):
1282 // ---------------------------------
1283 // [ 1 | 2 ]
1284 // [ --+-- ]
1285 // [ 3 | 4 ]
1286
1287 int i , j , k;
1288 int np1 = _nfree + 1;
1289 int nF = np1 + p1;
1290 double ** F = new double *[nF];
1291 double ** M = new double *[p1];
1292 for ( i = 0 ; i < nF ; ++i )
1293 F[i] = new double[nF];
1294
1295 // 1/4: MQ.MQ' (p+1,p+1):
1296 {
1297 for ( i = 0 ; i < p1 ; ++i )
1298 {
1299 M[i] = new double[_n_alpha];
1300 for ( j = 0 ; j < _n_alpha ; ++j )
1301 M[i][j] = compute_M ( i , j );
1302 for ( j = 0 ; j <= i ; ++j )
1303 {
1304 F[i][j] = 0.0;
1305 for ( k = np1 ; k < _n_alpha ; ++k )
1306 F[i][j] += M[i][k] * M[j][k];
1307 if ( i != j )
1308 F[j][i] = F[i][j];
1309 }
1310 }
1311 }
1312
1313 // 2/4: ML (p+1,n+1):
1314 for ( i = 0 ; i < p1 ; ++i )
1315 {
1316 F[i][p1] = 1.0;
1317 for ( j = p1+1 ; j < nF ; ++j )
1318 F[i][j] = M[i][j-p1];
1319 }
1320
1321 // 3/4: ML' (n+1,p+1):
1322 for ( j = 0 ; j < p1 ; ++j )
1323 {
1324 F[p1][j] = 1.0;
1325 for ( i = p1+1 ; i < nF ; ++i )
1326 F[i][j] = M[j][i-p1];
1327 }
1328
1329 // 4/4: 0 (n+1,n+1):
1330 for ( i = p1 ; i < nF ; ++i )
1331 for ( j = p1 ; j < nF ; ++j )
1332 F[i][j] = 0.0;
1333
1334
1335 #ifdef DEBUG
1336 _out << std::endl << "F=";
1337 for ( i = 0 ; i < nF ; ++i )
1338 {
1339 _out << "\t";
1340 for ( j = 0 ; j < nF ; ++j )
1341 _out << std::setw(12) << F[i][j] << " ";
1342 _out << std::endl;
1343 }
1344 #endif
1345
1346 for ( i = 0 ; i < p1 ; ++i )
1347 delete [] M[i];
1348 delete [] M;
1349
1350 bool error = false;
1351
1352 // SVD decomposition of the F matrix (F = U.W.V'):
1353 // -----------------------------------------------
1354 // (F will be transformed in U)
1355
1356 double * W = new double [nF];
1357 double ** V = new double *[nF];
1358 for ( i = 0 ; i < nF ; ++i )
1359 V[i] = new double[nF];
1360
1361 std::string error_msg;
1362
1363 if ( NOMAD::SVD_decomposition ( error_msg , F , W , V , nF , nF , max_mpn ) ) {
1364
1365 // compute condition number:
1366 compute_cond ( W , nF , eps );
1367
1368 #ifdef DEBUG
1369 _out << std::endl << "F=";
1370 for ( i = 0 ; i < nF ; ++i )
1371 {
1372 _out << "\t";
1373 for ( j = 0 ; j < nF ; ++j )
1374 _out << std::setw(12) << F[i][j] << " ";
1375 _out << std::endl;
1376 }
1377
1378 _out << std::endl << "W=\t";
1379 for ( i = 0 ; i < nF ; ++i )
1380 _out << std::setw(12) << W[i] << " ";
1381 _out << std::endl << std::endl << "cond=" << _cond << std::endl;
1382
1383 _out << std::endl << "V=";
1384 for ( i = 0 ; i < nF ; ++i )
1385 {
1386 _out << "\t";
1387 for ( j = 0 ; j < nF ; ++j )
1388 _out << std::setw(12) << V[i][j] << " ";
1389 _out << std::endl;
1390 }
1391 #endif
1392
1393 }
1394 else
1395 {
1396 #ifdef DEBUG
1397 _out << std::endl << "NOMAD::Quad_Model::construct_MFN_model(): "
1398 << "SVD decomposition (" << error_msg << ")"
1399 << std::endl << std::endl;
1400 #endif
1401 error = true;
1402 _cond.clear();
1403 }
1404
1405 // resolution of system F.[mu alpha_L]'=[f(Y) 0]' :
1406 // ------------------------------------------------
1407 if ( !error )
1408 {
1409 int m = static_cast<int> ( _bbot.size() );
1410 for ( i = 0 ; i < m ; ++i )
1411 if ( _alpha[i] )
1412 solve_MFN_system ( F , W , V , i , *_alpha[i] , eps );
1413 }
1414
1415 // free memory:
1416 for ( i = 0 ; i < nF ; ++i )
1417 {
1418 delete [] F[i];
1419 delete [] V[i];
1420 }
1421 delete [] F;
1422 delete [] V;
1423 delete [] W;
1424
1425 #ifdef DEBUG
1426 _out.close_block();
1427 #endif
1428
1429 return !error;
1430 }
1431
1432 /*--------------------------------------------------*/
1433 /* resolution of system F.[mu alpha_L]'=[f(Y) 0]' */
1434 /* for MFN interpolation (private) */
1435 /*--------------------------------------------------*/
solve_MFN_system(double ** F,double * W,double ** V,int bbo_index,NOMAD::Point & alpha,double eps) const1436 void NOMAD::Quad_Model::solve_MFN_system ( double ** F ,
1437 double * W ,
1438 double ** V ,
1439 int bbo_index ,
1440 NOMAD::Point & alpha ,
1441 double eps ) const
1442 {
1443 // resize the alpha vector:
1444 if ( alpha.size() != _n_alpha )
1445 alpha.reset ( _n_alpha , 0.0 );
1446
1447 int i , k , k1 , k2 ,
1448 np1 = _nfree + 1 ,
1449 nm1 = _nfree - 1 ,
1450 p1 = get_nY() ,
1451 nF = np1 + p1;
1452
1453 // step 1/2: find alpha_L and mu:
1454 // ---------
1455 double * alpha_tmp = new double [np1];
1456 double * mu_tmp = new double [ p1];
1457 double * mu = new double [ p1];
1458
1459 // if F is singular, some W values will be zero (or near zero);
1460 // each value that is smaller than eps is ignored:
1461 for ( i = 0 ; i < p1 ; ++i )
1462 {
1463 mu_tmp[i] = 0.0;
1464 if ( W[i] > eps )
1465 for ( k = 0 ; k < p1 ; ++k )
1466 mu_tmp[i] += F[k][i] *
1467 ( _Y[k]->get_bb_outputs()[bbo_index].value() ) / W[i];
1468 }
1469
1470 for ( i = p1 ; i < nF ; ++i )
1471 {
1472 alpha_tmp[i-p1] = 0.0;
1473 if ( W[i] > eps )
1474 for ( k = 0 ; k < p1 ; ++k )
1475 alpha_tmp[i-p1] += F[k][i] *
1476 ( _Y[k]->get_bb_outputs()[bbo_index].value() ) / W[i];
1477 }
1478
1479 for ( i = 0 ; i < p1 ; ++i )
1480 {
1481 mu[i] = 0.0;
1482 for ( k = 0 ; k < p1 ; ++k )
1483 mu[i] += V[i][k] * mu_tmp[k];
1484 for ( k = p1 ; k < nF ; ++k )
1485 mu[i] += V[i][k] * alpha_tmp[k-p1];
1486 }
1487
1488 for ( i = p1 ; i < nF ; ++i )
1489 {
1490 alpha[i-p1] = 0.0;
1491 for ( k = 0 ; k < p1 ; ++k )
1492 alpha[i-p1] += V[i][k] * mu_tmp[k];
1493 for ( k = p1 ; k < nF ; ++k )
1494 alpha[i-p1] += V[i][k] * alpha_tmp[k-p1];
1495 }
1496
1497 delete [] alpha_tmp;
1498 delete [] mu_tmp;
1499
1500 #ifdef DEBUG
1501 _out << std::endl << "output #" << bbo_index << ": mu=\t";
1502 for ( i = 0 ; i < p1 ; ++i )
1503 _out << std::setw(12) << mu[i] << " ";
1504 _out << std::endl;
1505
1506 _out << std::endl << "output #" << bbo_index << ": alpha_intermediate=\t";
1507 for ( i = 0 ; i < alpha.size() ; ++i )
1508 _out << std::setw(12) << alpha[i] << " ";
1509 _out << std::endl;
1510
1511 #endif
1512
1513 // step 2/2: find alpha_Q:
1514 // ---------
1515 for ( i = 0 ; i < _nfree ; ++i )
1516 {
1517 alpha[i+np1] = 0.0;
1518 for ( k = 0 ; k < p1 ; ++k )
1519 alpha[i+np1] += mu[k] * pow ( (*_Y[k])[_index[i+1]-1].value() , 2.0 ) / 2.0;
1520 }
1521
1522 for ( k1 = 0 ; k1 < nm1 ; ++k1 )
1523 for ( k2 = k1+1 ; k2 < _nfree ; ++k2 )
1524 {
1525 alpha[i+np1] = 0.0;
1526 for ( k = 0 ; k < p1 ; ++k )
1527 alpha[i+np1] += mu[k] * (*_Y[k])[_index[k1+1]-1].value() * (*_Y[k])[_index[k2+1]-1].value();
1528 ++i;
1529 }
1530
1531 delete [] mu;
1532 }
1533
1534 /*-----------------------------------------------------------*/
1535 /* check the interpolation set Y (private) */
1536 /*-----------------------------------------------------------*/
check_Y(void) const1537 bool NOMAD::Quad_Model::check_Y ( void ) const
1538 {
1539 if ( _Y.empty() ) {
1540 #ifdef DEBUG
1541 _out << std::endl << "NOMAD::Quad_Model::check_Y(): set Y is empty"
1542 << std::endl << std::endl;
1543 #endif
1544 return false;
1545 }
1546
1547 int nY = get_nY();
1548 int m = static_cast<int> ( _bbot.size() );
1549
1550 for ( int k = 0 ; k < nY ; ++k ) {
1551
1552 if ( _Y[k] == NULL ) {
1553 #ifdef DEBUG
1554 _out << std::endl
1555 << "NOMAD::Quad_Model::check_Y(): NULL pointer in the set Y"
1556 << std::endl << std::endl;
1557 #endif
1558 return false;
1559 }
1560
1561 if ( _Y[k]->get_eval_status() != NOMAD::EVAL_OK ) {
1562 #ifdef DEBUG
1563 _out << std::endl
1564 << "NOMAD::Quad_Model::check_Y(): a point in Y failed to evaluate"
1565 << std::endl << std::endl;
1566 #endif
1567 return false;
1568 }
1569
1570 const NOMAD::Point & bbo = _Y[k]->get_bb_outputs();
1571
1572 if ( !bbo.is_complete() ) {
1573 #ifdef DEBUG
1574 _out << std::endl
1575 << "NOMAD::Quad_Model::check_Y(): some bb outputs in Y are not defined"
1576 << std::endl << std::endl;
1577 #endif
1578 return false;
1579 }
1580
1581 if ( bbo.size() != m ) {
1582 #ifdef DEBUG
1583 _out << std::endl
1584 << "NOMAD::Quad_Model::check_Y(): "
1585 << "bb outputs in Y do not have the same dimension"
1586 << std::endl << std::endl;
1587 #endif
1588 return false;
1589 }
1590
1591 if ( _Y[k]->size() != _n ) {
1592 #ifdef DEBUG
1593 _out << std::endl
1594 << "NOMAD::Quad_Model::check_Y(): "
1595 << "points in Y do not have the same dimension"
1596 << std::endl << std::endl;
1597 #endif
1598 return false;
1599 }
1600 }
1601 return true;
1602 }
1603
1604 /*----------------------------------------------------*/
1605 /* check if the model is ready for evaluations */
1606 /*----------------------------------------------------*/
check(void) const1607 bool NOMAD::Quad_Model::check ( void ) const
1608 {
1609 if ( !_alpha )
1610 return false;
1611
1612 int nalpha = (_nfree+1)*(_nfree+2)/2;
1613 int i , m = static_cast<int> ( _bbot.size() );
1614
1615 for ( int bbo_index = 0 ; bbo_index < m ; ++bbo_index )
1616 {
1617
1618 if ( _alpha[bbo_index] )
1619 {
1620
1621 if ( _alpha[bbo_index]->size() != nalpha )
1622 return false;
1623
1624 for ( i = 0 ; i < nalpha ; ++i )
1625 if ( !(*_alpha[bbo_index])[i].is_defined() )
1626 return false;
1627 }
1628 }
1629
1630 return true;
1631 }
1632
1633
1634 /*--------------------------------------------------------------------------*/
1635 /* evaluate a model at a given point */
1636 /*--------------------------------------------------------------------------*/
1637 /* . the method assumes that x.size()==_n, alpha.is_complete(), and */
1638 /* alpha.size()==(_nfree+1)*(_nfree+2)/2 */
1639 /* . a more efficient version is used in Quad_Model_Evaluator::eval_x() ) */
1640 /*--------------------------------------------------------------------------*/
eval(const NOMAD::Point & x,const NOMAD::Point & alpha) const1641 NOMAD::Double NOMAD::Quad_Model::eval ( const NOMAD::Point & x ,
1642 const NOMAD::Point & alpha ) const
1643 {
1644 int i , j , k = 1 , nm1 = _n-1;
1645 NOMAD::Double z = alpha[0];
1646
1647 for ( i = 0 ; i < _n ; ++i )
1648 {
1649 if ( !_fixed_vars[i] )
1650 {
1651 z += x[i] * ( alpha[k] + 0.5 * alpha[k+_nfree] * x[i] );
1652 ++k;
1653 }
1654 }
1655
1656 k += _nfree;
1657
1658 for ( i = 0 ; i < nm1 ; ++i )
1659 if ( !_fixed_vars[i] )
1660 for ( j = i+1 ; j < _n ; ++j )
1661 if ( !_fixed_vars[j] )
1662 z += alpha[k++] * x[i] * x[j];
1663
1664 return z;
1665 }
1666
1667
1668 /*----------------------------------------------------------------*/
1669 /* compute model h and f values at a point */
1670 /*----------------------------------------------------------------*/
eval_hf(const NOMAD::Point & x,const NOMAD::Double & h_min,NOMAD::hnorm_type h_norm,NOMAD::Double & h,NOMAD::Double & f) const1671 void NOMAD::Quad_Model::eval_hf ( const NOMAD::Point & x ,
1672 const NOMAD::Double & h_min ,
1673 NOMAD::hnorm_type h_norm ,
1674 NOMAD::Double & h ,
1675 NOMAD::Double & f ) const
1676 {
1677 f.clear();
1678 h = 0.0;
1679
1680 int m = static_cast<int>(_bbot.size());
1681 NOMAD::Double bboi;
1682
1683 for ( int i = 0 ; i < m ; ++i ) {
1684
1685 if ( _alpha[i] ) {
1686
1687 bboi = eval ( x , *_alpha[i] );
1688
1689 if ( bboi.is_defined() ) {
1690
1691 if ( _bbot[i] == NOMAD::EB || _bbot[i] == NOMAD::PEB_E ) {
1692
1693 if ( bboi > h_min ) {
1694 h.clear();
1695 return;
1696 }
1697 }
1698
1699 else if ( ( _bbot[i] == NOMAD::FILTER ||
1700 _bbot[i] == NOMAD::PB ||
1701 _bbot[i] == NOMAD::PEB_P ) ) {
1702 if ( bboi > h_min ) {
1703 switch ( h_norm ) {
1704 case NOMAD::L1:
1705 h += bboi;
1706 break;
1707 case NOMAD::L2:
1708 h += bboi * bboi;
1709 break;
1710 case NOMAD::LINF:
1711 if ( bboi > h )
1712 h = bboi;
1713 break;
1714 }
1715 }
1716 }
1717
1718 else if ( _bbot[i] == NOMAD::OBJ )
1719 f = bboi;
1720 }
1721 }
1722 }
1723
1724 if ( h_norm == NOMAD::L2 )
1725 h = h.sqrt();
1726 }
1727
1728
1729
1730
1731 /*-----------------------------------------------------*/
1732 /* compute the maximal relative error of a model for */
1733 /* the interpolation set (private) */
1734 /*-----------------------------------------------------*/
compute_max_rel_err(void) const1735 NOMAD::Double NOMAD::Quad_Model::compute_max_rel_err ( void ) const
1736 {
1737 NOMAD::Double truth_value , model_value , rel_err , max_rel_err;
1738 int k , nY = get_nY() , m = static_cast<int> ( _bbot.size() );
1739
1740 for ( int bbo_index = 0 ; bbo_index < m ; ++bbo_index ) {
1741 if ( _alpha[bbo_index] ) {
1742 for ( k = 0 ; k < nY ; ++k ) {
1743 if ( _Y[k] && _Y[k]->get_eval_status() == NOMAD::EVAL_OK ) {
1744 truth_value = _Y[k]->get_bb_outputs()[bbo_index];
1745 if ( truth_value.is_defined() ) {
1746 model_value = eval ( *_Y[k] , *_alpha[bbo_index] );
1747 if ( model_value.is_defined() ) {
1748 if ( truth_value.abs() != 0.0 ) {
1749 rel_err = (truth_value-model_value).abs() / truth_value.abs();
1750 if ( !max_rel_err.is_defined() || rel_err > max_rel_err )
1751 max_rel_err = rel_err;
1752 }
1753 }
1754 }
1755 }
1756 }
1757 }
1758 }
1759 return max_rel_err;
1760 }
1761
1762 /*---------------------------------------------*/
1763 /* compute the cumulated error of a model */
1764 /* for the points of the interpolation set Y */
1765 /* and for one output in particular */
1766 /* (private) */
1767 /*---------------------------------------------*/
compute_model_error(int bbo_index,NOMAD::Double & error,NOMAD::Double & min_rel_err,NOMAD::Double & max_rel_err,NOMAD::Double & avg_rel_err) const1768 void NOMAD::Quad_Model::compute_model_error ( int bbo_index ,
1769 NOMAD::Double & error ,
1770 NOMAD::Double & min_rel_err ,
1771 NOMAD::Double & max_rel_err ,
1772 NOMAD::Double & avg_rel_err ) const
1773 {
1774 NOMAD::Double truth_value , model_value , rel_err;
1775 int nY = get_nY() , cnt = 0;
1776 bool chk = true;
1777
1778 max_rel_err.clear();
1779 min_rel_err.clear();
1780 avg_rel_err = error = 0.0;
1781
1782 #ifdef DEBUG
1783 std::ostringstream msg;
1784 msg << "output #" << bbo_index;
1785 _out.open_block ( msg.str() );
1786 #endif
1787
1788 for ( int k = 0 ; k < nY ; ++k )
1789 if ( _Y[k] && _Y[k]->get_eval_status() == NOMAD::EVAL_OK )
1790 {
1791 truth_value = _Y[k]->get_bb_outputs()[bbo_index];
1792
1793 if ( truth_value.is_defined() )
1794 {
1795 model_value = eval ( *_Y[k] , *_alpha[bbo_index] );
1796 if ( model_value.is_defined() )
1797 {
1798 rel_err.clear();
1799 if ( truth_value.abs() != 0.0 )
1800 rel_err = (truth_value-model_value).abs() / truth_value.abs();
1801 else
1802 {
1803 if (truth_value.abs()==model_value.abs())
1804 rel_err=0.0;
1805 else
1806 rel_err=NOMAD::INF;
1807 }
1808 if ( !max_rel_err.is_defined() || rel_err > max_rel_err )
1809 max_rel_err = rel_err;
1810 if ( !min_rel_err.is_defined() || rel_err < min_rel_err )
1811 min_rel_err = rel_err;
1812 avg_rel_err += rel_err;
1813 ++cnt;
1814
1815 #ifdef DEBUG
1816 _out << "Y[" << k << "]= ( ";
1817 _Y[k]->NOMAD::Point::display ( _out );
1818 _out << " )" << " f=" << truth_value
1819 << " m=" << model_value << " error^2="
1820 << ( model_value - truth_value ).pow2()
1821 << " rel_err=" << rel_err
1822 << std::endl;
1823 #endif
1824 error += ( model_value - truth_value ).pow2();
1825 }
1826 else
1827 {
1828 chk = false;
1829 break;
1830 }
1831 }
1832 else
1833 {
1834 chk = false;
1835 break;
1836 }
1837 }
1838
1839 #ifdef DEBUG
1840 _out.close_block();
1841 #endif
1842
1843 if ( chk)
1844 { // Case where chk is true (at least one model_value and the corresponding thruth value were defined => cnt != 0)
1845 error = error.sqrt();
1846 avg_rel_err = avg_rel_err / cnt;
1847 }
1848 else
1849 {
1850 error.clear();
1851 min_rel_err.clear();
1852 max_rel_err.clear();
1853 avg_rel_err.clear();
1854 }
1855 }
1856
1857 /*-----------------------------------------------------------*/
1858 /* display the model coefficients */
1859 /*-----------------------------------------------------------*/
display_model_coeffs(const NOMAD::Display & out) const1860 void NOMAD::Quad_Model::display_model_coeffs ( const NOMAD::Display & out ) const
1861 {
1862 if ( _error_flag ) {
1863 out << "model coefficients: could not be constructed" << std::endl;
1864 return;
1865 }
1866
1867 int m = static_cast<int> ( _bbot.size() );
1868
1869 out << NOMAD::open_block ( "model coefficients" );
1870 for ( int i = 0 ; i < m ; ++i ) {
1871 out << "output #";
1872 out.display_int_w ( i , m );
1873 out << ": ";
1874 if ( _alpha[i] ) {
1875 out<< "[ ";
1876 _alpha[i]->display ( out , " " , 6 );
1877 out << " ]";
1878 }
1879 else
1880 out << "NULL";
1881 out << std::endl;
1882 }
1883 out.close_block();
1884 }
1885
1886 /*-----------------------------------------------------------*/
1887 /* display the interpolation set Y */
1888 /*-----------------------------------------------------------*/
display_Y(const NOMAD::Display & out,const std::string & title) const1889 void NOMAD::Quad_Model::display_Y ( const NOMAD::Display & out ,
1890 const std::string & title ) const
1891 {
1892 out << NOMAD::open_block ( title );
1893 int nY = get_nY();
1894 for ( int k = 0 ; k < nY ; ++k ) {
1895 out << "#";
1896 out.display_int_w ( k , nY );
1897 out << ": ";
1898 if ( _Y[k] ) {
1899 out << "( ";
1900 _Y[k]->NOMAD::Point::display ( out , " " , 12 );
1901 out << " ) bbo=[ ";
1902 _Y[k]->get_bb_outputs().display ( out , " " , 12 );
1903 out << " ]";
1904 }
1905 else
1906 out << "NULL";
1907 out << std::endl;
1908 }
1909 out.close_block();
1910 }
1911
1912 /*-------------------------------------------------------*/
1913 /* display cumulated error on the interpolation points */
1914 /*-------------------------------------------------------*/
display_Y_error(const NOMAD::Display & out) const1915 void NOMAD::Quad_Model::display_Y_error ( const NOMAD::Display & out ) const
1916 {
1917 if ( _error_flag ) {
1918 out << "model error on the interpolation set: cannot be computed"
1919 << std::endl;
1920 return;
1921 }
1922
1923 int i ;
1924 int index = -1;
1925 int m = static_cast<int> ( _bbot.size() );
1926
1927 for ( i = 0 ; i < m ; ++i )
1928 if ( _alpha[i] ) {
1929 if ( index >= 0 ) {
1930 index = -1;
1931 break;
1932 }
1933 else
1934 index = i;
1935 }
1936
1937 NOMAD::Double error , min_rel_err , max_rel_err , avg_rel_err;
1938
1939 // only one output:
1940 if ( index >= 0 ) {
1941 compute_model_error ( index , error , min_rel_err , max_rel_err , avg_rel_err );
1942 out << "model errors on the interpolation set: error="
1943 << error << " min_rel_err=" << min_rel_err
1944 << " max_rel_err=" << max_rel_err << " avg_rel_err=" << avg_rel_err
1945 << std::endl;
1946 }
1947
1948 // several outputs:
1949 else {
1950
1951 out.open_block ( "model error on the interpolation set" );
1952
1953 NOMAD::Double error_i , min_rel_err_i , max_rel_err_i , avg_rel_err_i;
1954
1955 error = avg_rel_err = 0.0;
1956 min_rel_err.clear();
1957 max_rel_err.clear();
1958
1959 int cnt = 0;
1960
1961 for ( i = 0 ; i < m ; ++i )
1962 if ( _alpha[i] ) {
1963
1964 ++cnt;
1965
1966 compute_model_error ( i ,
1967 error_i ,
1968 min_rel_err_i ,
1969 max_rel_err_i ,
1970 avg_rel_err_i );
1971
1972 if (error_i.is_defined())
1973 error += error_i;
1974 if (avg_rel_err_i.is_defined())
1975 avg_rel_err += avg_rel_err_i;
1976 if ( !min_rel_err.is_defined() || min_rel_err_i < min_rel_err )
1977 min_rel_err = min_rel_err_i;
1978 if ( !max_rel_err.is_defined() || max_rel_err_i > max_rel_err )
1979 max_rel_err = max_rel_err_i;
1980
1981 out << "output #" << i << ": error=" << error_i
1982 << " min_rel_err=" << min_rel_err_i << " max_rel_err="
1983 << max_rel_err_i << " avg_rel_err=" << avg_rel_err_i << std::endl;
1984 }
1985
1986 out << std::endl << "global: error=" << error
1987 << " min_rel_err=" << min_rel_err
1988 << " max_rel_err=" << max_rel_err
1989 << " avg_rel_err=" << avg_rel_err / cnt
1990 << std::endl << NOMAD::close_block();
1991 }
1992 }
1993
1994 /*-----------------------------------------------*/
1995 /* display Lagrange polynomials (private) */
1996 /*-----------------------------------------------*/
display_lagrange_polynomials(const std::vector<NOMAD::Point * > & l,const std::vector<NOMAD::Eval_Point * > & Y) const1997 void NOMAD::Quad_Model::display_lagrange_polynomials
1998 ( const std::vector<NOMAD::Point *> & l ,
1999 const std::vector<NOMAD::Eval_Point *> & Y ) const
2000 {
2001 int i , j , nY = static_cast<int> ( Y.size() );
2002
2003 // display Lagrange polynomials:
2004 _out << std::endl << NOMAD::open_block ( "Lagrange polynomials" );
2005 for ( i = 0 ; i < _n_alpha ; ++i ) {
2006 _out << "l[";
2007 _out.display_int_w ( i , _n_alpha );
2008 _out << "] = [ ";
2009 l[i]->NOMAD::Point::display ( _out , " " , 14 , -1 );
2010 _out << "]" << std::endl;
2011 }
2012 _out.close_block();
2013
2014 // display current set Y:
2015 _out << std::endl << NOMAD::open_block ( "current set Y" );
2016 for ( i = 0 ; i < nY ; ++i ) {
2017 _out << "Y[";
2018 _out.display_int_w ( i , nY );
2019 _out << "] = ";
2020 if ( Y[i] ) {
2021 _out << "( ";
2022 Y[i]->NOMAD::Point::display ( _out , " " , 6 , -1 );
2023 _out << " )";
2024 }
2025 else
2026 _out << "NULL";
2027 _out << std::endl;
2028 }
2029
2030 _out.close_block();
2031
2032 // display l(Y): should be the identity matrix:
2033 NOMAD::Double tmp , err = 0.0;
2034 _out << std::endl << NOMAD::open_block ( "l(Y)" );
2035 for ( i = 0 ; i < _n_alpha ; ++i ) {
2036 _out << "l[";
2037 _out.display_int_w ( i , _n_alpha );
2038 _out << "]: ";
2039 for ( j = 0 ; j < _n_alpha ; ++j ) {
2040 tmp.clear();
2041 if ( j < nY && Y[j] ) {
2042 tmp = eval ( *Y[j] , *l[i] );
2043 err += (i==j) ? (tmp-1.0).abs() : tmp.abs();
2044 }
2045 tmp.display ( _out , "%6.3f" );
2046 _out << " ";
2047 }
2048 _out << std::endl;
2049 }
2050 _out << std::endl << "error (with identity) = "
2051 << err << std::endl << NOMAD::close_block() << std::endl;
2052 }
2053