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 Directions.cpp
38 \brief Set of polling directions (implementation)
39 \author Sebastien Le Digabel
40 \date 2010-04-13
41 \see Directions.hpp
42 */
43 #include "Directions.hpp"
44 #include <math.h>
45 #include "XMesh.hpp"
46 #include "SMesh.hpp"
47
48
49 /*---------------------------------------------------------*/
50 /* constructor */
51 /*---------------------------------------------------------*/
Directions(int nc,const std::set<NOMAD::direction_type> & direction_types,const std::set<NOMAD::direction_type> & sec_poll_dir_types,const NOMAD::Display & out)52 NOMAD::Directions::Directions
53 ( int nc ,
54 const std::set<NOMAD::direction_type> & direction_types ,
55 const std::set<NOMAD::direction_type> & sec_poll_dir_types ,
56 const NOMAD::Display & out )
57 : _nc ( nc ) ,
58 _direction_types ( direction_types ) ,
59 _sec_poll_dir_types ( sec_poll_dir_types ) ,
60 _is_binary ( false ) ,
61 _is_categorical ( false ) ,
62 _lt_initialized ( false ) ,
63 _out ( out )
64 {
65 // check the directions:
66 if ( _direction_types.find ( NOMAD::NO_DIRECTION ) != _direction_types.end() )
67 _direction_types.clear();
68
69 if ( _sec_poll_dir_types.find ( NOMAD::NO_DIRECTION ) != _sec_poll_dir_types.end() )
70 _sec_poll_dir_types.clear();
71
72 // is_orthomads: true if at least one direction is of type Ortho-MADS:
73 _is_orthomads = NOMAD::dirs_have_orthomads ( _direction_types );
74 if ( !_is_orthomads )
75 _is_orthomads = NOMAD::dirs_have_orthomads ( _sec_poll_dir_types );
76
77
78 }
79
80 /*---------------------------------------------------------*/
81 /* destructor */
82 /*---------------------------------------------------------*/
~Directions(void)83 NOMAD::Directions::~Directions ( void )
84 {
85 if ( _lt_initialized )
86 {
87 int n = 2 * NOMAD::L_LIMITS;
88 for ( int i = 0 ; i <= n ; ++i )
89 delete _bl[i];
90 }
91
92 }
93
94 /*---------------------------------------------------------*/
95 /* LT-MADS initializations (private) */
96 /*---------------------------------------------------------*/
lt_mads_init(void)97 void NOMAD::Directions::lt_mads_init ( void )
98 {
99 int n = 2 * NOMAD::L_LIMITS;
100 for ( int i = 0 ; i <= n ; ++i )
101 {
102 _bl [i] = NULL;
103 _hat_i[i] = -1;
104 }
105 _lt_initialized = true;
106 }
107
108 /*---------------------------------------------------------*/
109 /* set_binary */
110 /*---------------------------------------------------------*/
set_binary(void)111 void NOMAD::Directions::set_binary ( void )
112 {
113 _is_binary = true;
114 _is_categorical = false;
115 _is_orthomads = false;
116 _direction_types.clear();
117 _direction_types.insert ( NOMAD::GPS_BINARY );
118 if ( !_sec_poll_dir_types.empty() ) {
119 _sec_poll_dir_types.clear();
120 _sec_poll_dir_types.insert ( NOMAD::GPS_BINARY );
121 }
122 }
123
124 /*---------------------------------------------------------*/
125 /* set_categorical */
126 /*---------------------------------------------------------*/
set_categorical(void)127 void NOMAD::Directions::set_categorical ( void )
128 {
129 _is_categorical = true;
130 _is_binary = false;
131 _is_orthomads = false;
132 _direction_types.clear();
133 _sec_poll_dir_types.clear();
134 }
135
136 /*----------------------------------------------------------------------*/
137 /* compute binary directions when all groups of variables are binary */
138 /* (private) */
139 /*----------------------------------------------------------------------*/
compute_binary_directions(std::list<NOMAD::Direction> & d) const140 void NOMAD::Directions::compute_binary_directions
141 ( std::list<NOMAD::Direction> & d ) const
142 {
143 // _GPS_BINARY_ n directions:
144 NOMAD::Direction * pd;
145 for ( int i = 0 ; i < _nc ; ++i ) {
146 d.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_BINARY ) );
147 pd = &(*(--d.end()));
148 (*pd)[i] = 1.0;
149 }
150 }
151
152
153 /*----------------------------------------------------------------------------*/
154 /* get the directions on a unit n-sphere for a given mesh */
155 /*----------------------------------------------------------------------------*/
compute(std::list<NOMAD::Direction> & dirs,NOMAD::poll_type poll,const NOMAD::OrthogonalMesh & mesh)156 void NOMAD::Directions::compute ( std::list<NOMAD::Direction> & dirs ,
157 NOMAD::poll_type poll ,
158 const NOMAD::OrthogonalMesh & mesh )
159 {
160
161 dirs.clear();
162
163 // categorical variables: do nothing:
164 if ( _is_categorical )
165 return;
166
167 // binary variables: we use special directions:
168 if ( _is_binary )
169 {
170 compute_binary_directions ( dirs );
171 return;
172 }
173
174 NOMAD::Double pow_gps , cst;
175 const NOMAD::Direction * bl;
176 NOMAD::Direction * pd;
177 int i , j , k , hat_i;
178 std::list<NOMAD::Direction>::iterator it2 , end2;
179
180 const NOMAD::Point mesh_indices=mesh.get_mesh_indices();
181 int mesh_index=int(mesh_indices[0].value()); // single mesh_index used only for isotropic mesh
182
183
184 /*-----------------------------------*/
185 /* loop on the types of directions */
186 /*-----------------------------------*/
187 const std::set<NOMAD::direction_type> & dir_types =
188 (poll == NOMAD::PRIMARY) ? _direction_types : _sec_poll_dir_types;
189
190 std::set<NOMAD::direction_type>::const_iterator it , end = dir_types.end() ;
191 for ( it = dir_types.begin() ; it != end ; ++it )
192 {
193
194 if ( *it == NOMAD::UNDEFINED_DIRECTION ||
195 *it == NOMAD::NO_DIRECTION ||
196 *it == NOMAD::MODEL_SEARCH_DIR )
197 continue;
198
199 /*--------------------------------------------------*/
200 /* Ortho-MADS */
201 /* ---> dirs are on a unit n-sphere by construction */
202 /*--------------------------------------------------*/
203 if ( NOMAD::dir_is_orthomads (*it) )
204 {
205 bool success_dir;
206
207
208 NOMAD::Direction dir ( _nc , 0.0 , *it );
209 NOMAD::Double alpha_t_l;
210
211 success_dir=compute_dir_on_unit_sphere ( dir );
212
213 if ( success_dir )
214 {
215 #ifdef DEBUG
216 if ( dynamic_cast<const NOMAD::XMesh*> ( &mesh ) )
217 {
218 _out << std::endl
219 << NOMAD::open_block ( "compute Ortho-MADS directions with XMesh" )
220 << "type = " << *it << std::endl;
221 NOMAD::Point mesh_indices=mesh.get_mesh_indices();
222 _out<< "Mesh indices = ( ";
223 mesh_indices.NOMAD::Point::display( _out, " " , -1 , -1 );
224 _out << " )" << std::endl;
225 _out<< "Unit sphere Direction = ( ";
226 dir.NOMAD::Point::display ( _out , " " , -1 , -1 );
227 _out << " )" << std::endl << NOMAD::close_block();
228 }
229 else if ( dynamic_cast<const NOMAD::SMesh*> ( &mesh ) )
230 {
231 _out << std::endl
232 << NOMAD::open_block ( "compute Ortho-MADS directions with SMesh" )
233 << "type = " << *it
234 << " on isotropic mesh " << std::endl
235 << "mesh index (lk) = " << mesh_index << std::endl
236 << "alpha (tk,lk) = " << alpha_t_l << std::endl
237 << "Unit sphere Direction = ( ";
238 dir.NOMAD::Point::display ( _out , " " , -1 , -1 );
239 _out << " )" << std::endl << NOMAD::close_block();
240
241 }
242
243 #endif
244 // Ortho-MADS 2n and n+1:
245 // ----------------------
246 if ( *it == NOMAD::ORTHO_2N || *it == NOMAD::ORTHO_NP1_QUAD || *it == NOMAD::ORTHO_NP1_NEG )
247 {
248
249 // creation of the 2n directions:
250 // For ORTHO_NP1 the reduction from 2N to N+1 is done in MADS::POLL
251 NOMAD::Direction ** H = new NOMAD::Direction * [2*_nc];
252
253 // Ordering D_k alternates Hk and -Hk instead of [H_k -H_k]
254 for ( i = 0 ; i < _nc ; ++i )
255 {
256 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , *it ) );
257 H[i ] = &(*(--dirs.end()));
258 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , *it ) );
259 H[i+_nc] = &(*(--dirs.end()));
260 }
261
262 // Householder transformations on the 2n directions on a unit n-sphere:
263 householder ( dir , true , H );
264
265
266 delete [] H;
267
268
269 }
270
271 // Ortho-MADS 1 or Ortho-MADS 2:
272 // -----------------------------
273 else
274 {
275 dirs.push_back ( dir );
276 if ( *it == NOMAD::ORTHO_2 )
277 dirs.push_back ( -dir );
278 }
279 }
280 }
281
282 /*--------------------------------*/
283 /* LT-MADS */
284 /* ---> dirs NOT on unit n-sphere */
285 /*--------------------------------*/
286 else if ( NOMAD::dir_is_ltmads (*it) )
287 {
288
289 if ( !_lt_initialized)
290 lt_mads_init();
291
292 bl = get_bl ( mesh , *it , hat_i );
293
294 // LT-MADS 1 or LT-MADS 2: -b(l) and/or b(l):
295 // ------------------------------------------
296 if ( *it == NOMAD::LT_1 || *it == NOMAD::LT_2 )
297 {
298 dirs.push_back ( - *bl );
299 if ( *it == NOMAD::LT_2 )
300 dirs.push_back ( *bl );
301 }
302
303 // LT-MADS 2n or LT-MADS n+1:
304 // --------------------------
305 else {
306
307 // row permutation vector:
308 int * row_permutation_vector = new int [_nc];
309 row_permutation_vector[hat_i] = hat_i;
310
311 NOMAD::Random_Pickup rp ( _nc );
312
313 for ( i = 0 ; i < _nc ; ++i )
314 if ( i != hat_i )
315 {
316 j = rp.pickup();
317 if ( j == hat_i )
318 j = rp.pickup();
319 row_permutation_vector[i] = j;
320 }
321
322 rp.reset();
323
324 for ( j = 0 ; j < _nc ; ++j )
325 {
326 i = rp.pickup();
327 if ( i != hat_i )
328 {
329
330 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , *it ) );
331 pd = &(*(--dirs.end()));
332
333 create_lt_direction ( mesh , *it , i , hat_i , pd );
334 permute_coords ( *pd , row_permutation_vector );
335 }
336 else
337 dirs.push_back ( *bl );
338
339 // completion to maximal basis:
340 if ( *it == NOMAD::LT_2N )
341 dirs.push_back ( NOMAD::Direction ( - *(--dirs.end()) , NOMAD::LT_2N ) );
342
343 }
344
345 delete [] row_permutation_vector;
346
347 // completion to minimal basis:
348 if ( *it == NOMAD::LT_NP1 )
349 {
350 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::LT_NP1 ) );
351 pd = &(*(--dirs.end()));
352
353 end2 = --dirs.end();
354 for ( it2 = dirs.begin() ; it2 != end2 ; ++it2 )
355 {
356 for ( i = 0 ; i < _nc ; ++i )
357 (*pd)[i] -= (*it2)[i];
358 }
359 }
360 }
361
362
363 }
364
365 /*--------------------------------*/
366 /* GPS */
367 /* ---> dirs NOT on unit n-sphere */
368 /*--------------------------------*/
369 else
370 {
371
372 // GPS for binary variables: should'nt be here:
373 if ( *it == NOMAD::GPS_BINARY )
374 {
375 #ifdef DEBUG
376 _out << std::endl << "Warning (" << "Directions.cpp" << ", " << __LINE__
377 << "): tried to generate binary directions at the wrong place)"
378 << std::endl << std::endl;
379 #endif
380 dirs.clear();
381 return;
382 }
383
384 // this value represents the non-zero values in GPS directions
385 // (it is tau^{|ell_k|/2}, and it is such that delta_k * pow_gps = Delta_k):
386 if ( !pow_gps.is_defined() )
387 pow_gps = pow ( mesh.get_update_basis() , abs(mesh_index) / 2.0 );
388
389 // GPS 2n, static:
390 // ---------------
391
392 if ( *it == NOMAD::GPS_2N_STATIC )
393 {
394 for ( i = 0 ; i < _nc ; ++i )
395 {
396 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_2N_STATIC ) );
397 pd = &(*(--dirs.end()));
398 (*pd)[i] = pow_gps;
399
400 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_2N_STATIC ) );
401 pd = &(*(--dirs.end()));
402 (*pd)[i] = -pow_gps;
403 }
404 }
405
406 // GPS 2n, random:
407 // ---------------
408 else if ( *it == NOMAD::GPS_2N_RAND )
409 {
410
411 int v , cnt;
412
413 std::list <int>::const_iterator end3;
414 std::list <int>::iterator it3;
415 std::list <int> rem_cols;
416 std::vector<int> rem_col_by_row ( _nc );
417
418 // creation of the 2n directions:
419 std::vector<NOMAD::Direction *> pdirs ( 2 * _nc );
420
421 for ( i = 0 ; i < _nc ; ++i )
422 {
423
424 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_2N_RAND ) );
425 pdirs[i] = &(*(--dirs.end()));
426 (*pdirs[i])[i] = pow_gps;
427
428 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_2N_RAND ) );
429 pdirs[i+_nc] = &(*(--dirs.end()));
430 (*pdirs[i+_nc])[i] = -pow_gps;
431
432 rem_cols.push_back(i);
433 rem_col_by_row[i] = i;
434 }
435
436 // random perturbations:
437 for ( i = 1 ; i < _nc ; ++i )
438 {
439 if ( rem_col_by_row[i] > 0 )
440 {
441 v = NOMAD::RNG::rand()%3 - 1; // v in { -1 , 0 , 1 }
442 if ( v )
443 {
444
445 // we decide a (i,j) couple:
446 k = NOMAD::RNG::rand()%(rem_col_by_row[i]);
447 cnt = 0;
448 end3 = rem_cols.end();
449 it3 = rem_cols.begin();
450 while ( cnt != k )
451 {
452 ++it3;
453 ++cnt;
454 }
455 j = *it3;
456
457 // the perturbation:
458 (*pdirs[i])[j] = (*pdirs[j+_nc])[i] = -v * pow_gps;
459 (*pdirs[j])[i] = (*pdirs[i+_nc])[j] = v * pow_gps;
460
461 // updates:
462 rem_cols.erase(it3);
463 it3 = rem_cols.begin();
464 while ( *it3 != i )
465 ++it3;
466 rem_cols.erase(it3);
467 for ( k = i+1 ; k < _nc ; ++k )
468 rem_col_by_row[k] -= j<k ? 2 : 1;
469 }
470 }
471 }
472 }
473
474 // GPS n+1, static:
475 // ----------------
476 else if ( *it == NOMAD::GPS_NP1_STATIC )
477 {
478
479 // (n+1)^st direction:
480 dirs.push_back ( NOMAD::Direction ( _nc , -pow_gps , NOMAD::GPS_NP1_STATIC ) );
481
482 // first n directions:
483 for ( i = 0 ; i < _nc ; ++i )
484 {
485 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_STATIC ) );
486 pd = &(*(--dirs.end()));
487 (*pd)[i] = pow_gps;
488 }
489 }
490
491 // GPS n+1, random:
492 // ----------------
493 else if ( *it == NOMAD::GPS_NP1_RAND )
494 {
495
496 // (n+1)^st direction:
497 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_RAND ) );
498 NOMAD::Direction * pd1 = &(*(--dirs.end()));
499
500 NOMAD::Double d;
501
502 // first n directions:
503 for ( i = 0 ; i < _nc ; ++i )
504 {
505 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_RAND ) );
506 pd = &(*(--dirs.end()));
507
508 d = NOMAD::RNG::rand()%2 ? -pow_gps : pow_gps;
509 (*pd )[i] = d;
510 (*pd1)[i] = -d;
511 }
512 }
513
514 // GPS n+1, static, uniform angles:
515 // --------------------------------
516 else if ( *it == NOMAD::GPS_NP1_STATIC_UNIFORM )
517 {
518
519 cst = pow_gps * sqrt(static_cast<double>(_nc)*(_nc+1))/_nc;
520
521 // n first directions:
522 for ( j = _nc-1 ; j >= 0 ; --j )
523 {
524 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_STATIC_UNIFORM ) );
525 pd = &(*(--dirs.end()));
526
527 k = _nc-j-1;
528
529 for ( i = 0 ; i < k ; ++i )
530 (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));
531
532 (*pd)[k] = (cst * (j+1)) / sqrt(static_cast<double>(j+1)*(j+2));
533
534 }
535
536 // (n+1)^st direction:
537 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_STATIC_UNIFORM ) );
538 pd = &(*(--dirs.end()));
539 for ( i = 0 ; i < _nc ; ++i )
540 (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));
541
542 }
543
544 // GPS n+1, random, uniform angles:
545 // --------------------------------
546 // (we apply the procedure defined in
547 // "Pattern Search Methods for user-provided points:
548 // application to molecular geometry problems",
549 // by Alberto, Nogueira, Rocha and Vicente,
550 // SIOPT 14-4, 1216-1236, 2004, doi:10.1137/S1052623400377955)
551 else if ( *it == NOMAD::GPS_NP1_RAND_UNIFORM )
552 {
553
554 cst = pow_gps * sqrt(static_cast<double>(_nc)*(_nc+1))/_nc;
555
556 // n first directions:
557 for ( j = _nc-1 ; j >= 0 ; --j )
558 {
559 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_STATIC_UNIFORM ) );
560 pd = &(*(--dirs.end()));
561
562 k = _nc-j-1;
563
564 for ( i = 0 ; i < k ; ++i )
565 (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));
566
567 (*pd)[k] = (cst * (j+1)) / sqrt(static_cast<double>(j+1)*(j+2));
568
569 }
570
571 // (n+1)^st direction:
572 dirs.push_back ( NOMAD::Direction ( _nc , 0.0 , NOMAD::GPS_NP1_STATIC_UNIFORM ) );
573 pd = &(*(--dirs.end()));
574 for ( i = 0 ; i < _nc ; ++i )
575 (*pd)[i] = -cst / sqrt(static_cast<double>(_nc-i)*(_nc-i+1));
576
577 // random perturbations:
578 NOMAD::Random_Pickup rp ( _nc );
579 std::vector<bool> done ( _nc );
580 bool chg_sgn;
581 std::list<NOMAD::Direction>::iterator it;
582 NOMAD::Double tmp;
583
584 end2 = dirs.end();
585 for ( i = 0 ; i < _nc ; ++i )
586 done[i] = false;
587
588 for ( i = 0 ; i < _nc ; ++i )
589 {
590 k = rp.pickup();
591 if ( i != k && !done[i] )
592 {
593 chg_sgn = ( NOMAD::RNG::rand()%2 != 0 );
594
595 for ( it = dirs.begin() ; it != end2 ; ++it )
596 {
597 tmp = (*it)[i];
598 (*it)[i] = ( chg_sgn ? -1.0 : 1.0 ) * (*it)[k];
599 (*it)[k] = tmp;
600 }
601
602 done[i] = done[k] = true;
603 }
604 else
605 if ( NOMAD::RNG::rand()%2 )
606 for ( it = dirs.begin() ; it != end2 ; ++it )
607 (*it)[i] = -(*it)[i];
608 }
609 }
610 }
611 } // end of loop on the types of directions
612
613 // The direction are unscaled on a unit n-sphere when necessary (that is directions that are not orthomads)
614 for ( it2 = dirs.begin() ; it2 != dirs.end() ; ++it2 )
615 if ( ! NOMAD::dir_is_orthomads ( it2->get_type() ) )
616 {
617 NOMAD::Double norm_dir=it2->norm();
618 for (int i=0 ; i < _nc ; ++i )
619 (*it2)[i] /= norm_dir;
620 }
621
622 }
623
624
625 /*-----------------------------------------------------------------------------*/
626 /* compute a random direction on a unit N-Sphere */
627 /* see http://en.wikipedia.org/wiki/N-sphere */
628 /*-----------------------------------------------------------------------------*/
compute_dir_on_unit_sphere(NOMAD::Direction & random_dir) const629 bool NOMAD::Directions::compute_dir_on_unit_sphere ( NOMAD::Direction & random_dir ) const
630 {
631
632 int i;
633 NOMAD::Double norm,d;
634
635 for ( i = 0 ; i < _nc ; ++i )
636 random_dir[i]=NOMAD::RNG::normal_rand(0,1);
637
638 norm=random_dir.norm();
639
640 if ( norm==0 )
641 return false;
642
643
644 for ( i = 0 ; i < _nc ; ++i )
645 random_dir[i]/=norm;
646
647 return true;
648
649 }
650
651
652 /*-----------------------------------------------------------------*/
653 /* compute the squared norm of normalized(2u_t-e) for Ortho-MADS */
654 /* (private) */
655 /*-----------------------------------------------------------------*/
eval_ortho_norm(const NOMAD::Double & x,const NOMAD::Double & norm,const NOMAD::Point & b,NOMAD::Point & new_b) const656 NOMAD::Double NOMAD::Directions::eval_ortho_norm ( const NOMAD::Double & x ,
657 const NOMAD::Double & norm ,
658 const NOMAD::Point & b ,
659 NOMAD::Point & new_b ) const
660 {
661 NOMAD::Double fx = 0.0;
662
663 for ( int i = 0 ; i < _nc ; ++i ) {
664 new_b[i] = ( x * b[i] / norm ).round();
665 fx += new_b[i]*new_b[i];
666 }
667
668 return fx;
669 }
670
671 /*--------------------------------------------------------*/
672 /* get the expression of an integer t in inverse base p */
673 /* (private, static) */
674 /*--------------------------------------------------------*/
get_phi(int t,int p)675 NOMAD::Double NOMAD::Directions::get_phi ( int t , int p )
676 {
677 int div;
678 int size = int ( ceil ( log(static_cast<double>(t+1)) /
679 log(static_cast<double>(p)) ) );
680 int ll = t;
681 NOMAD::Double d = 0.0;
682
683 for ( int i = 0 ; i < size ; ++i ) {
684 div = NOMAD::Double ( pow ( p , size-i-1.0 ) ).round();
685 d += ( ll / div ) * pow ( static_cast<double>(p) , i-size );
686 ll = ll % div;
687 }
688 return d;
689 }
690
691 /*----------------------------------------------------------------*/
692 /* . Householder transformation to generate _nc directions from */
693 /* a given direction */
694 /* . compute also H[i+nc] = -H[i] (completion to 2n directions) */
695 /* . private method */
696 /*----------------------------------------------------------------*/
householder(const NOMAD::Direction & dir,bool complete_to_2n,NOMAD::Direction ** H) const697 void NOMAD::Directions::householder ( const NOMAD::Direction & dir ,
698 bool complete_to_2n ,
699 NOMAD::Direction ** H ) const
700 {
701 int i , j;
702
703 NOMAD::Double norm2 = dir.squared_norm() , v , h2i;
704
705 for ( i = 0 ; i < _nc ; ++i ) {
706
707 h2i = 2 * dir[i];
708
709 for ( j = 0 ; j < _nc ; ++j ) {
710
711 // H[i]:
712 (*H[i])[j] = v = (i==j) ? norm2 - h2i * dir[j] : - h2i * dir[j];
713
714 // -H[i]:
715 if ( complete_to_2n )
716 (*H[i+_nc])[j] = -v;
717 }
718 }
719 }
720
721 /*---------------------------------------------------------*/
722 /* get the LT-MADS b(l) direction (private) */
723 /*---------------------------------------------------------*/
get_bl(const NOMAD::OrthogonalMesh & mesh,NOMAD::direction_type dtype,int & hat_i)724 const NOMAD::Direction * NOMAD::Directions::get_bl ( const NOMAD::OrthogonalMesh & mesh,
725 NOMAD::direction_type dtype ,
726 int & hat_i )
727 {
728
729 NOMAD::Point mesh_indices = mesh.get_mesh_indices();
730 int mesh_index =static_cast<int>(mesh_indices[0].value());
731 NOMAD::Direction * bl = _bl [ mesh_index + NOMAD::L_LIMITS ];
732 hat_i = _hat_i [ mesh_index + NOMAD::L_LIMITS ];
733
734 if ( !bl )
735 {
736 hat_i = -1;
737 create_lt_direction ( mesh , dtype , -1 , hat_i , bl );
738 }
739
740 return bl;
741 }
742
743 /*---------------------------------------------------------*/
744 /* create a new LT-MADS direction (private) */
745 /*---------------------------------------------------------*/
746 /* (if hat_i == -1, a new b(l) direction */
747 /* is created and hat_i is set) */
748 /*---------------------------------------------------------*/
create_lt_direction(const NOMAD::OrthogonalMesh & mesh,NOMAD::direction_type dtype,int diag_i,int & hat_i,NOMAD::Direction * & dir)749 void NOMAD::Directions::create_lt_direction ( const NOMAD::OrthogonalMesh & mesh ,
750 NOMAD::direction_type dtype ,
751 int diag_i ,
752 int & hat_i ,
753 NOMAD::Direction *& dir )
754 {
755
756 int mesh_index=static_cast<int>(mesh.get_mesh_indices()[0].value());
757 int i_pow_tau =
758 static_cast<int>
759 ( ceil ( pow ( mesh.get_update_basis() , abs(mesh_index) / 2.0 ) ) );
760
761 int j = diag_i+1;
762
763 // new b(l) direction:
764 if ( hat_i < 0 )
765 {
766 _hat_i [ mesh_index + NOMAD::L_LIMITS ] = diag_i = hat_i = NOMAD::RNG::rand()%_nc;
767 _bl [ mesh_index + NOMAD::L_LIMITS ] = dir
768 = new NOMAD::Direction ( _nc, 0.0, dtype );
769
770 j = 0;
771 }
772
773 (*dir)[diag_i] = (NOMAD::RNG::rand()%2) ? -i_pow_tau : i_pow_tau;
774
775 for ( int k = j ; k < _nc ; ++k )
776 if ( k != hat_i )
777 {
778 (*dir)[k] = NOMAD::RNG::rand()%i_pow_tau;
779 if ( NOMAD::RNG::rand()%2 && (*dir)[k] > 0.0 )
780 (*dir)[k] = -(*dir)[k];
781 }
782
783 #ifdef DEBUG
784 if ( j==0 )
785 _out << "new LT-MADS b(l) direction: b(" << mesh_index << ") = "
786 << *dir << std::endl << std::endl;
787 #endif
788 }
789
790 /*---------------------------------------------------------*/
791 /* permute the coordinates of a direction (private) */
792 /*---------------------------------------------------------*/
permute_coords(NOMAD::Direction & dir,const int * permutation_vector) const793 void NOMAD::Directions::permute_coords ( NOMAD::Direction & dir ,
794 const int * permutation_vector ) const
795 {
796 NOMAD::Point tmp = dir;
797 for ( int i = 0 ; i < _nc ; ++i )
798 dir [ permutation_vector[i] ] = tmp[i];
799 }
800
801
802 /*---------------------------------------------------------*/
803 /* display */
804 /*---------------------------------------------------------*/
display(const NOMAD::Display & out) const805 void NOMAD::Directions::display ( const NOMAD::Display & out ) const
806 {
807 out << "n : " << _nc << std::endl
808 << "types : { ";
809 std::set<NOMAD::direction_type>::const_iterator it , end = _direction_types.end();
810 for ( it = _direction_types.begin() ; it != end ; ++it )
811 out << "[" << *it << "] ";
812 out << "}" << std::endl
813 << "sec poll types: { ";
814 end = _sec_poll_dir_types.end();
815 for ( it = _sec_poll_dir_types.begin() ; it != end ; ++it )
816 out << "[" << *it << "] ";
817 out << "}" << std::endl;
818 if ( _is_orthomads )
819 {
820 out << "seed : " << NOMAD::RNG::get_seed() << endl;
821 }
822 }
823
824 /*---------------------------------------------------------*/
825 /* comparison operator */
826 /*---------------------------------------------------------*/
operator <(const NOMAD::Directions & d) const827 bool NOMAD::Directions::operator < ( const NOMAD::Directions & d ) const
828 {
829 // number of variables:
830 if ( _nc < d._nc )
831 return true;
832 if ( d._nc < _nc )
833 return false;
834
835 // boolean variables:
836 if ( _is_binary != d._is_binary )
837 return _is_binary;
838 if ( _is_categorical != d._is_categorical )
839 return _is_categorical;
840 if ( _is_orthomads != d._is_orthomads )
841 return _is_orthomads;
842
843 // direction types:
844 size_t nd = _direction_types.size();
845 if ( nd < d._direction_types.size() )
846 return true;
847 if ( d._direction_types.size() < nd )
848 return false;
849
850 size_t ns = _sec_poll_dir_types.size();
851 if ( ns < d._sec_poll_dir_types.size() )
852 return true;
853 if ( d._sec_poll_dir_types.size() < ns )
854 return false;
855
856 std::set<NOMAD::direction_type>::const_iterator
857 it1 = _direction_types.begin() ,
858 it2 = d._direction_types.begin() ,
859 end = _direction_types.end();
860
861 while ( it1 != end ) {
862 if ( *it1 < *it2 )
863 return true;
864 if ( *it2 < *it1 )
865 return false;
866 ++it1;
867 ++it2;
868 }
869
870 it1 = _sec_poll_dir_types.begin();
871 it2 = d._sec_poll_dir_types.begin();
872 end = _sec_poll_dir_types.end();
873
874 while ( it1 != end ) {
875 if ( *it1 < *it2 )
876 return true;
877 if ( *it2 < *it1 )
878 return false;
879 ++it1;
880 ++it2;
881 }
882
883 return false;
884 }
885
886