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 utils.cpp
38 \brief Utility functions
39 \author Sebastien Le Digabel
40 \date 2010-03-23
41 \see utils.hpp
42 */
43 #include "utils.hpp"
44
45 /*---------------------------------------------------------------*/
46 /* construct the first n prime numbers */
47 /*---------------------------------------------------------------*/
construct_primes(int n,int * primes)48 void NOMAD::construct_primes ( int n , int * primes )
49 {
50 bool is_prime;
51 double kk;
52 int i = 0 , k = 2 , j;
53 while ( true ) {
54 is_prime = true;
55 kk = sqrt ( static_cast<double>(k) );
56 for ( j = 2 ; j <= kk ; ++j )
57 if ( k%j==0 ) {
58 is_prime = false;
59 break;
60 }
61 if ( is_prime ) {
62 primes[i++] = k;
63 if ( i==n )
64 break;
65 }
66 ++k;
67 }
68 }
69
70 /*--------------------------------------------*/
71 /* decompose a string (sentence) into words */
72 /*--------------------------------------------*/
get_words(const std::string & sentence,std::list<std::string> & words)73 void NOMAD::get_words ( const std::string & sentence , std::list<std::string> & words )
74 {
75 std::string s;
76 std::istringstream in ( sentence );
77 while ( true ) {
78 in >> s;
79 if ( in.fail() )
80 break;
81 words.push_back ( s );
82 }
83 }
84
85 /*---------------------------------------------------------------*/
86 /* get the pid (process id): used as random seed or unique tag */
87 /*---------------------------------------------------------------*/
get_pid(void)88 int NOMAD::get_pid ( void )
89 {
90 #ifdef _MSC_VER
91 return _getpid();
92 #else
93 return getpid();
94 #endif
95 }
96
97 /*---------------------------------------*/
98 /* called at the beginning of NOMAD */
99 /*---------------------------------------*/
begin(int argc,char ** argv)100 void NOMAD::begin ( int argc , char ** argv )
101 {
102 #ifdef USE_MPI
103 MPI_Init ( &argc , &argv );
104 #endif
105 }
106
107 /*---------------------------------------*/
108 /* called at the end of NOMAD */
109 /*---------------------------------------*/
end(void)110 void NOMAD::end ( void )
111 {
112 #ifdef USE_MPI
113 MPI_Finalize();
114 #endif
115 }
116
117 /*-----------------------------------------------------------------*/
118 /* check if a file exists and is executable */
119 /*-----------------------------------------------------------------*/
check_exe_file(const std::string & file_name)120 bool NOMAD::check_exe_file ( const std::string & file_name )
121 {
122 #ifdef WINDOWS
123 // don't check on Windows:
124 return true;
125 #else
126 return ( access ( file_name.c_str() , X_OK ) == 0 );
127 #endif
128 }
129
130 /*-----------------------------------------------------------------*/
131 /* check if a file exists and is readable */
132 /*-----------------------------------------------------------------*/
check_read_file(const std::string & file_name)133 bool NOMAD::check_read_file ( const std::string & file_name )
134 {
135 #ifdef _MSC_VER
136 return ( _access ( file_name.c_str() , 4 ) == 0 );
137 #else
138 return ( access ( file_name.c_str() , R_OK ) == 0 );
139 #endif
140 }
141
142 /*-----------------------------------------------------------------*/
143 /* NOMAD::itos */
144 /*-----------------------------------------------------------------*/
itos(int i)145 std::string NOMAD::itos ( int i )
146 {
147 std::ostringstream oss;
148 oss << i;
149 return oss.str();
150 }
151
152 /*-----------------------------------------------------------------*/
153 /* NOMAD::itos */
154 /*-----------------------------------------------------------------*/
itos(size_t i)155 std::string NOMAD::itos ( size_t i )
156 {
157 std::ostringstream oss;
158 oss << i;
159 return oss.str();
160 }
161
162
163 /*-----------------------------------------------------------------*/
164 /* NOMAD::toupper - 1/2 */
165 /*-----------------------------------------------------------------*/
toupper(std::string & s)166 void NOMAD::toupper ( std::string & s )
167 {
168 size_t ns = s.size();
169 for ( size_t i = 0 ; i < ns ; ++i )
170 s[i] = std::toupper(s[i]);
171 }
172
173 /*-----------------------------------------------------------------*/
174 /* NOMAD::toupper - 2/2 */
175 /*-----------------------------------------------------------------*/
toupper(std::list<std::string> & ls)176 void NOMAD::toupper ( std::list<std::string> & ls )
177 {
178 std::list<std::string>::iterator it;
179 std::list<std::string>::const_iterator end = ls.end();
180 for ( it = ls.begin() ; it != end ; ++it )
181 NOMAD::toupper ( *it );
182 }
183
184 /*-----------------------------------------------------------------*/
185 /* NOMAD::atoi */
186 /*-----------------------------------------------------------------*/
atoi(const std::string & s,int & i)187 bool NOMAD::atoi ( const std::string & s , int & i )
188 {
189 i = -1;
190 if ( s.empty() )
191 return false;
192
193 size_t n = s.size();
194
195 if ( s[0] == '-' ) {
196 if ( n > 1 && s[1] == '-' )
197 return false;
198 std::string ss = s;
199 ss.erase(ss.begin());
200 if ( NOMAD::atoi ( ss , i ) ) {
201 i = -i;
202 return true;
203 }
204 return false;
205 }
206
207 for ( size_t k = 0 ; k < n ; ++k )
208 if ( !isdigit(s[k]) )
209 return false;
210 i = std::atoi(s.c_str());
211 return true;
212 }
213
atoi(char c,int & i)214 bool NOMAD::atoi ( char c , int & i ) {
215 std::string s = "-";
216 s[0] = c;
217 return NOMAD::atoi(s,i);
218 }
219
220 /*-------------------------------------------------------------------*/
221 /* if a NOMAD::bb_output_type variable corresponds to a constraint */
222 /*-------------------------------------------------------------------*/
bbot_is_constraint(NOMAD::bb_output_type bbot)223 bool NOMAD::bbot_is_constraint ( NOMAD::bb_output_type bbot )
224 {
225 return ( bbot == NOMAD::EB ||
226 bbot == NOMAD::PB ||
227 bbot == NOMAD::PEB_P ||
228 bbot == NOMAD::PEB_E ||
229 bbot == NOMAD::FILTER );
230 }
231
232 /*-----------------------------------------------------------------------*/
233 /* if a NOMAD::direction_type variable corresponds to a MADS direction */
234 /*-----------------------------------------------------------------------*/
dir_is_mads(NOMAD::direction_type dt)235 bool NOMAD::dir_is_mads ( NOMAD::direction_type dt )
236 {
237 return (dt == NOMAD::ORTHO_1 ||
238 dt == NOMAD::ORTHO_2 ||
239 dt == NOMAD::ORTHO_NP1_QUAD ||
240 dt == NOMAD::ORTHO_NP1_NEG ||
241 dt == NOMAD::DYN_ADDED ||
242 dt == NOMAD::ORTHO_2N ||
243 dt == NOMAD::LT_1 ||
244 dt == NOMAD::LT_2 ||
245 dt == NOMAD::LT_2N ||
246 dt == NOMAD::LT_NP1 );
247 }
248
249 /*----------------------------------------------------------------------*/
250 /* if a NOMAD::direction_type variable corresponds to a GPS direction */
251 /*----------------------------------------------------------------------*/
dir_is_gps(NOMAD::direction_type dt)252 bool NOMAD::dir_is_gps ( NOMAD::direction_type dt )
253 {
254 return ( dt == NOMAD::GPS_BINARY ||
255 dt == NOMAD::GPS_2N_STATIC ||
256 dt == NOMAD::GPS_2N_RAND ||
257 dt == NOMAD::GPS_NP1_STATIC_UNIFORM ||
258 dt == NOMAD::GPS_NP1_STATIC ||
259 dt == NOMAD::GPS_NP1_RAND_UNIFORM ||
260 dt == NOMAD::GPS_NP1_RAND );
261 }
262
263 /*--------------------------------------------------------------------------*/
264 /* if a NOMAD::direction_type variable corresponds to a LT-MADS direction */
265 /*--------------------------------------------------------------------------*/
dir_is_ltmads(NOMAD::direction_type dt)266 bool NOMAD::dir_is_ltmads ( NOMAD::direction_type dt )
267 {
268 return ( dt == NOMAD::LT_1 ||
269 dt == NOMAD::LT_2 ||
270 dt == NOMAD::LT_2N ||
271 dt == NOMAD::LT_NP1 );
272 }
273
274 /*-------------------------------------------------------------------------*/
275 /* if a NOMAD::direction_type variable corresponds to a random direction */
276 /*-------------------------------------------------------------------------*/
dir_is_random(NOMAD::direction_type dt)277 bool NOMAD::dir_is_random ( NOMAD::direction_type dt )
278 {
279 return ( dt == NOMAD::GPS_NP1_RAND ||
280 dt == NOMAD::GPS_NP1_RAND_UNIFORM ||
281 dt == NOMAD::GPS_2N_RAND ||
282 dt == NOMAD::LT_1 ||
283 dt == NOMAD::LT_2 ||
284 dt == NOMAD::LT_2N ||
285 dt == NOMAD::LT_NP1 );
286 }
287
288 /*-----------------------------------------------------------------------------*/
289 /* if a NOMAD::direction_type variable corresponds to a Ortho-MADS direction */
290 /*-----------------------------------------------------------------------------*/
dir_is_orthomads(NOMAD::direction_type dt)291 bool NOMAD::dir_is_orthomads ( NOMAD::direction_type dt )
292 {
293 return (dt == NOMAD::ORTHO_1 ||
294 dt == NOMAD::ORTHO_2 ||
295 dt == NOMAD::ORTHO_NP1_QUAD ||
296 dt == NOMAD::ORTHO_NP1_NEG ||
297 dt == NOMAD::ORTHO_2N );
298 }
299
300 /*---------------------------------------------------------------------*/
301 /* check if a set of directions includes one Ortho-MADS direction */
302 /* (true if at least one direction in the set is of type Ortho-MADS) */
303 /*---------------------------------------------------------------------*/
dirs_have_orthomads(const std::set<NOMAD::direction_type> & dir_types)304 bool NOMAD::dirs_have_orthomads ( const std::set<NOMAD::direction_type> & dir_types )
305 {
306 std::set<NOMAD::direction_type>::const_iterator it , end = dir_types.end();
307 for ( it = dir_types.begin() ; it != end ; ++it )
308 if ( NOMAD::dir_is_orthomads (*it) )
309 return true;
310 return false;
311 }
312
313 /*---------------------------------------------------------------------*/
314 /* check if a set of directions includes one Ortho-MADS N+1 direction */
315 /* (true if at least one direction in the set is of type Ortho-MADS) */
316 /*---------------------------------------------------------------------*/
dirs_have_orthomads_np1(const std::set<NOMAD::direction_type> & dir_types)317 bool NOMAD::dirs_have_orthomads_np1( const std::set<NOMAD::direction_type> & dir_types )
318 {
319 std::set<NOMAD::direction_type>::const_iterator it , end = dir_types.end();
320 for ( it = dir_types.begin() ; it != end ; ++it )
321 if ( (*it)==NOMAD::ORTHO_NP1_QUAD || (*it)==NOMAD::ORTHO_NP1_NEG)
322 return true;
323 return false;
324 }
325
326 /*---------------------------------------------------*/
327 /* returns true if one of the string of ls is in s */
328 /*---------------------------------------------------*/
string_find(const std::string & s,const std::list<std::string> & ls)329 bool NOMAD::string_find ( const std::string & s , const std::list<std::string> & ls )
330 {
331 std::list<std::string>::const_iterator it , end = ls.end();
332 for ( it = ls.begin() ; it != end ; ++it )
333 if ( NOMAD::string_find ( s , *it ) )
334 return true;
335 return false;
336 }
337
338 /*---------------------------------------------------*/
339 /* search a string into another string */
340 /*---------------------------------------------------*/
string_find(const std::string & s1,const std::string & s2)341 bool NOMAD::string_find ( const std::string & s1 , const std::string & s2 )
342 {
343 return ( s1.find(s2) < s1.size() );
344 }
345
346 /*--------------------------------------------------------------------------*/
347 /* returns true if one of the string in ls is an element of s */
348 /*--------------------------------------------------------------------------*/
string_match(const std::string & s,const std::list<std::string> & ls)349 bool NOMAD::string_match ( const std::string & s , const std::list<std::string> & ls )
350 {
351 std::list<std::string>::const_iterator it , end = ls.end();
352 for ( it = ls.begin() ; it != end ; ++it )
353 if ( s.compare(*it) == 0 )
354 return true;
355 return false;
356 }
357
358
359 /*-----------------------------------------------------------------*/
360 /* interpret a list of strings as a direction type */
361 /*-----------------------------------------------------------------*/
strings_to_direction_type(const std::list<std::string> & ls,NOMAD::direction_type & dt)362 bool NOMAD::strings_to_direction_type ( const std::list<std::string> & ls ,
363 NOMAD::direction_type & dt )
364 {
365
366 dt = NOMAD::UNDEFINED_DIRECTION;
367
368 if ( ls.empty() || ls.size() > 4 )
369 return false;
370
371 std::list<std::string>::const_iterator it = ls.begin() , end = ls.end();
372 std::string s = *it;
373 NOMAD::toupper ( s );
374
375 // no direction:
376 if ( s == "NONE" ) {
377 dt = NOMAD::NO_DIRECTION;
378 return true;
379 }
380
381 // Ortho-MADS with 1, 2, n+1 (plus QUAD or NEG), or 2n directions:
382 if ( s == "ORTHO" )
383 {
384 ++it;
385 if ( it == end )
386 {
387 dt = NOMAD::ORTHO_NP1_QUAD; // Default for ORTHO
388 return true;
389 }
390 if ( *it == "1" )
391 {
392 dt = NOMAD::ORTHO_1;
393 return true;
394 }
395 if ( *it == "2" )
396 {
397 dt = NOMAD::ORTHO_2;
398 return true;
399 }
400 s = *it;
401 NOMAD::toupper ( s );
402 if ( s == "2N" )
403 {
404 dt = NOMAD::ORTHO_2N;
405 return true;
406 }
407 if ( s == "N+1" )
408 {
409 ++it;
410 if (it==end)
411 {
412 dt = NOMAD::ORTHO_NP1_QUAD; // Default for ORTHO N+1
413 return true;
414 }
415 s = *it;
416 NOMAD::toupper ( s );
417 if ( s=="QUAD" )
418 {
419 dt= NOMAD::ORTHO_NP1_QUAD;
420 return true;
421 }
422 if ( s=="NEG" )
423 {
424 dt=NOMAD::ORTHO_NP1_NEG;
425 return true;
426 }
427
428 }
429
430 return false;
431 }
432
433 // LT-MADS with 1, 2 or 2n directions:
434 if ( s == "LT" ) {
435 ++it;
436 if ( it == end ) {
437 dt = NOMAD::LT_2N;
438 return true;
439 }
440 if ( *it == "1" ) {
441 dt = NOMAD::LT_1;
442 return true;
443 }
444 if ( *it == "2" ) {
445 dt = NOMAD::LT_2;
446 return true;
447 }
448 s = *it;
449 NOMAD::toupper ( s );
450 if ( s == "N+1" ) {
451 dt = NOMAD::LT_NP1;
452 return true;
453 }
454 if ( s == "2N" ) {
455 dt = NOMAD::LT_2N;
456 return true;
457 }
458 return false;
459 }
460
461 // GPS:
462 if ( s == "GPS" ) {
463 ++it;
464 if ( it == end ) {
465 dt = NOMAD::GPS_2N_STATIC;
466 return true;
467 }
468 s = *it;
469 NOMAD::toupper ( s );
470
471 // GPS for binary variables:
472 if ( s == "BINARY" || s == "BIN" ) {
473 dt = NOMAD::GPS_BINARY;
474 return true;
475 }
476
477 // GPS, n+1 directions:
478 if ( s == "N+1" ) {
479 ++it;
480 if ( it == end ) {
481 dt = NOMAD::GPS_NP1_STATIC;
482 return true;
483 }
484 s = *it;
485 NOMAD::toupper ( s );
486
487 // GPS, n+1, static:
488 if ( s == "STATIC" ) {
489 ++it;
490 if ( it == end ) {
491 dt = NOMAD::GPS_NP1_STATIC;
492 return true;
493 }
494 s = *it;
495 NOMAD::toupper ( s );
496 if ( s == "UNIFORM" ) {
497 dt = NOMAD::GPS_NP1_STATIC_UNIFORM;
498 return true;
499 }
500 return false;
501 }
502
503 // GPS, n+1, random:
504 if ( s == "RAND" || s == "RANDOM" ) {
505 ++it;
506 if ( it == end ) {
507 dt = NOMAD::GPS_NP1_RAND;
508 return true;
509 }
510 s = *it;
511 NOMAD::toupper ( s );
512 if ( s == "UNIFORM" ) {
513 dt = NOMAD::GPS_NP1_RAND_UNIFORM;
514 return true;
515 }
516 return false;
517 }
518 return false;
519 }
520
521 // 2n directions:
522 if ( s == "2N" ) {
523 ++it;
524 if ( it == end ) {
525 dt = NOMAD::GPS_2N_STATIC;
526 return true;
527 }
528 s = *it;
529 NOMAD::toupper ( s );
530 if ( s == "STATIC" ) {
531 dt = NOMAD::GPS_2N_STATIC;
532 return true;
533 }
534 if ( s == "RAND" || s == "RANDOM" ) {
535 dt = NOMAD::GPS_2N_RAND;
536 return true;
537 }
538 return false;
539 }
540 return false;
541 }
542 return false;
543 }
544
545 /*---------------------------------------*/
546 /* convert a string into a hnorm_type */
547 /*---------------------------------------*/
string_to_hnorm_type(const std::string & s,NOMAD::hnorm_type & hn)548 bool NOMAD::string_to_hnorm_type ( const std::string & s , NOMAD::hnorm_type & hn )
549 {
550 std::string ss = s;
551 NOMAD::toupper(ss);
552 if ( ss == "L1" ) {
553 hn = NOMAD::L1;
554 return true;
555 }
556 if ( ss == "L2" ) {
557 hn = NOMAD::L2;
558 return true;
559 }
560 if ( ss == "LINF" ) {
561 hn = NOMAD::LINF;
562 return true;
563 }
564 return false;
565 }
566
567 /*-----------------------------------------*/
568 /* convert a string into a TGP_mode_type */
569 /*-----------------------------------------*/
string_to_TGP_mode_type(const std::string & s,NOMAD::TGP_mode_type & m)570 bool NOMAD::string_to_TGP_mode_type ( const std::string & s , NOMAD::TGP_mode_type & m )
571 {
572 std::string ss = s;
573 NOMAD::toupper(ss);
574 if ( ss == "FAST" ) {
575 m = NOMAD::TGP_FAST;
576 return true;
577 }
578 if ( ss == "PRECISE" ) {
579 m = NOMAD::TGP_PRECISE;
580 return true;
581 }
582 if ( ss == "USER" ) {
583 m = NOMAD::TGP_USER;
584 return true;
585 }
586 return false;
587 }
588
589 /*--------------------------------------------------*/
590 /* convert a string into a multi_formulation_type */
591 /*--------------------------------------------------*/
string_to_multi_formulation_type(const std::string & s,NOMAD::multi_formulation_type & mft)592 bool NOMAD::string_to_multi_formulation_type ( const std::string & s ,
593 NOMAD::multi_formulation_type & mft )
594 {
595 std::string ss = s;
596 NOMAD::toupper(ss);
597 if ( ss == "NORMALIZED" ) {
598 mft = NOMAD::NORMALIZED;
599 return true;
600 }
601 if ( ss == "PRODUCT" ) {
602 mft = NOMAD::PRODUCT;
603 return true;
604 }
605 if ( ss == "DIST_L1" ) {
606 mft = NOMAD::DIST_L1;
607 return true;
608 }
609 if ( ss == "DIST_L2" ) {
610 mft = NOMAD::DIST_L2;
611 return true;
612 }
613 if ( ss == "DIST_LINF" ) {
614 mft = NOMAD::DIST_LINF;
615 return true;
616 }
617 return false;
618 }
619
620 /*-------------------------------------------*/
621 /* convert a string into a bb_output_type */
622 /*-------------------------------------------*/
string_to_bb_output_type(const std::string & s,NOMAD::bb_output_type & bbot)623 bool NOMAD::string_to_bb_output_type ( const std::string & s ,
624 NOMAD::bb_output_type & bbot )
625 {
626 std::string ss = s;
627 NOMAD::toupper(ss);
628
629 if ( ss == "OBJ" ) {
630 bbot = NOMAD::OBJ;
631 return true;
632 }
633 if ( ss == "EB" ) {
634 bbot = NOMAD::EB;
635 return true;
636 }
637 if ( ss == "PB" || ss == "CSTR" ) {
638 bbot = NOMAD::PB;
639 return true;
640 }
641 if ( ss == "PEB" ) {
642 bbot = NOMAD::PEB_P;
643 return true;
644 }
645 if ( ss == "F" ) {
646 bbot = NOMAD::FILTER;
647 return true;
648 }
649 if ( ss == "STAT_AVG" ) {
650 bbot = NOMAD::STAT_AVG;
651 return true;
652 }
653 if ( ss == "STAT_SUM" ) {
654 bbot = NOMAD::STAT_SUM;
655 return true;
656 }
657 if ( ss == "CNT_EVAL" ) {
658 bbot = NOMAD::CNT_EVAL;
659 return true;
660 }
661 if ( ss == "NOTHING" || ss == "-" ) {
662 bbot = NOMAD::UNDEFINED_BBO;
663 return true;
664 }
665 return false;
666 }
667
668 /*-----------------------------------------------------------------*/
669 /* convert a string into a bb_input_type */
670 /*-----------------------------------------------------------------*/
string_to_bb_input_type(const std::string & s,NOMAD::bb_input_type & bbit)671 bool NOMAD::string_to_bb_input_type ( const std::string & s ,
672 NOMAD::bb_input_type & bbit )
673 {
674 std::string ss = s;
675 NOMAD::toupper ( ss );
676 if ( ss=="R" || ss=="REAL" ) {
677 bbit = NOMAD::CONTINUOUS;
678 return true;
679 }
680 if ( ss=="C" || ss=="CAT" ) {
681 bbit = NOMAD::CATEGORICAL;
682 return true;
683 }
684 if ( ss=="B" || ss=="BIN" ) {
685 bbit = NOMAD::BINARY;
686 return true;
687 }
688 if ( ss=="I" || ss=="INT" ) {
689 bbit = NOMAD::INTEGER;
690 return true;
691 }
692 return false;
693 }
694
695 /*-----------------------------------------------------------------*/
696 /* convert a string into a model_type */
697 /*-----------------------------------------------------------------*/
string_to_model_type(const std::string & s,NOMAD::model_type & mt)698 bool NOMAD::string_to_model_type ( const std::string & s ,
699 NOMAD::model_type & mt )
700 {
701 std::string ss = s;
702 NOMAD::toupper ( ss );
703 if ( ss=="TGP" || ss=="TGP_MODEL" ) {
704 mt = NOMAD::TGP_MODEL;
705 return true;
706 }
707 if ( ss=="QUADRATIC" || ss=="QUADRATIC_MODEL" ) {
708 mt = NOMAD::QUADRATIC_MODEL;
709 return true;
710 }
711
712 mt = NOMAD::NO_MODEL;
713 return false;
714 }
715
716 /*----------------------------------------------------------------------*/
717 /* convert a string in {"YES","NO","Y","N"} to a bool */
718 /* value of return: -1: error */
719 /* 0: bool=false */
720 /* 1: bool=true */
721 /*----------------------------------------------------------------------*/
string_to_bool(const std::string & ss)722 int NOMAD::string_to_bool ( const std::string & ss )
723 {
724 std::string s = ss;
725 NOMAD::toupper ( s );
726 if ( s=="Y" || s=="YES" || s=="1" || s=="TRUE" )
727 return 1;
728 if ( s=="N" || s=="NO" || s=="0" || s=="FALSE" )
729 return 0;
730 return -1;
731 }
732
733 /*---------------------------------------------------*/
734 /* convert a string 'i-j' to the integers i and j */
735 /*---------------------------------------------------*/
string_to_index_range(const std::string & s,int & i,int & j,int * n,bool check_order)736 bool NOMAD::string_to_index_range ( const std::string & s ,
737 int & i ,
738 int & j ,
739 int * n ,
740 bool check_order )
741 {
742 if ( s.empty() )
743 return false;
744
745 if ( s == "*" ) {
746 if ( !n )
747 return false;
748 i = 0;
749 j = *n-1;
750 return true;
751 }
752
753 if ( s[0] == '-' ) {
754
755 size_t ns = s.size();
756 if ( ns > 1 && s[1] == '-' )
757 return false;
758
759 std::string ss = s;
760 ss.erase ( ss.begin() );
761
762 if ( NOMAD::string_to_index_range ( ss , i , j , n , false ) ) {
763 i = -i;
764 return true;
765 }
766 return false;
767 }
768
769 std::istringstream in (s);
770 std::string s1;
771
772 getline ( in , s1 , '-' );
773
774 if (in.fail())
775 return false;
776
777 size_t k , n1 = s1.size();
778
779 if ( n1 >= s.size() - 1 ) {
780 for ( k = 0 ; k < n1 ; ++k )
781 if (!isdigit(s1[k]))
782 return false;
783 if ( !NOMAD::atoi ( s1 , i ) )
784 return false;
785 if ( n1 == s.size() ) {
786 j = i;
787 return true;
788 }
789 if (n) {
790 j = *n-1;
791 return true;
792 }
793 return false;
794 }
795
796 std::string s2;
797 getline (in, s2);
798
799 if (in.fail())
800 return false;
801
802 size_t n2 = s2.size();
803 for ( k = 0 ; k < n2 ; ++k )
804 if ( !isdigit(s2[k]) )
805 return false;
806
807 if ( !NOMAD::atoi ( s1, i ) || !NOMAD::atoi ( s2 , j ) )
808 return false;
809
810 return !check_order || i <= j;
811 }
812
813
get_rank(double ** M,size_t m,size_t n)814 int NOMAD::get_rank(double ** M,
815 size_t m,
816 size_t n)
817 {
818 double * W = new double [n];
819 double ** V = new double *[n];
820 for (size_t i = 0 ; i < n ; ++i )
821 {
822 V[i]=new double [n];
823 }
824
825 std::string error_msg;
826 NOMAD::SVD_decomposition ( error_msg , M , W , V , static_cast<int>(m) , static_cast<int>(n) );
827
828 for (size_t i=0;i<n;++i)
829 delete [] V[i];
830 delete [] V;
831
832
833 if (! error_msg.empty())
834 {
835 delete [] W;
836 return -1;
837 }
838
839 int rank=0;
840 for (size_t i=0;i<n;i++)
841 {
842 if (fabs(W[i])>NOMAD::SVD_EPS)
843 rank++;
844 }
845
846 delete [] W;
847 return rank;
848
849 }
850
851
852 /*--------------------------------------------------------------*/
853 /* SVD decomposition */
854 /* inspired and recoded from an old numerical recipes version */
855 /*--------------------------------------------------------------*/
856 /* */
857 /* M = U . W . V' */
858 /* */
859 /* M ( m x n ) */
860 /* U ( m x n ) */
861 /* W ( n x n ) */
862 /* V ( n x n ) */
863 /* */
864 /* U.U' = U'.U = I if m = n */
865 /* U'.U = I if m > m */
866 /* */
867 /* V.V' = V'.V = I */
868 /* */
869 /* W diagonal */
870 /* */
871 /* M is given as first argument and becomes U */
872 /* W is given as a size-n vector */
873 /* V is given, not V' */
874 /* */
875 /*--------------------------------------------------------------*/
SVD_decomposition(std::string & error_msg,double ** M,double * W,double ** V,int m,int n,int max_mpn)876 bool NOMAD::SVD_decomposition ( std::string & error_msg ,
877 double ** M ,
878 double * W ,
879 double ** V ,
880 int m ,
881 int n ,
882 int max_mpn ) // default=1500
883 {
884 error_msg.clear();
885
886 if ( max_mpn > 0 && m+n > max_mpn ) {
887 error_msg = "SVD_decomposition() error: m+n > " + NOMAD::itos ( max_mpn );
888 return false;
889 }
890
891 double * rv1 = new double[n];
892 double scale = 0.0;
893 double g = 0.0;
894 double norm = 0.0;
895
896 int nm1 = n - 1;
897
898 bool flag;
899 int i , j , k , l , its , jj , nm = 0;
900 double s , f , h , tmp , c , x , y , z , absf , absg , absh;
901
902 const int NITER = 30;
903
904 // Initialization W and V
905 for (i=0; i < n; ++i)
906 {
907 W[i]=0.0;
908 for (j=0; j < n ; ++j)
909 V[i][j]=0.0;
910 }
911
912 // Householder reduction to bidiagonal form:
913 for ( i = 0 ; i < n ; ++i )
914 {
915 l = i + 1;
916 rv1[i] = scale * g;
917 g = s = scale = 0.0;
918 if ( i < m )
919 {
920 for ( k = i ; k < m ; ++k )
921 scale += fabs ( M[k][i] );
922 if ( scale ) {
923 for ( k = i ; k < m ; ++k )
924 {
925 M[k][i] /= scale;
926 s += M[k][i] * M[k][i];
927 }
928 f = M[i][i];
929 g = ( f >= 0.0 ) ? -fabs(sqrt(s)) : fabs(sqrt(s));
930 h = f * g - s;
931 M[i][i] = f - g;
932 for ( j = l ; j < n ; ++j ) {
933 for ( s = 0.0 , k = i ; k < m ; ++k )
934 s += M[k][i] * M[k][j];
935 f = s / h;
936 for ( k = i ; k < m ; ++k )
937 M[k][j] += f * M[k][i];
938 }
939 for ( k = i ; k < m ; ++k )
940 M[k][i] *= scale;
941 }
942 }
943 W[i] = scale * g;
944 g = s = scale = 0.0;
945 if ( i < m && i != nm1 )
946 {
947 for ( k = l ; k < n ; ++k )
948 scale += fabs ( M[i][k] );
949 if ( scale )
950 {
951 for ( k = l ; k < n ; ++k )
952 {
953 M[i][k] /= scale;
954 s += M[i][k] * M[i][k];
955 }
956 f = M[i][l];
957 g = ( f >= 0.0 ) ? -fabs(sqrt(s)) : fabs(sqrt(s));
958 h = f * g - s;
959 M[i][l] = f - g;
960 for ( k = l ; k < n ; ++k )
961 rv1[k] = M[i][k] / h;
962 for ( j = l ; j < m ; ++j ) {
963 for ( s=0.0,k=l ; k < n ; ++k )
964 s += M[j][k] * M[i][k];
965 for ( k=l ; k < n ; ++k )
966 M[j][k] += s * rv1[k];
967 }
968 for ( k = l ; k < n ; ++k )
969 M[i][k] *= scale;
970 }
971 }
972 tmp = fabs ( W[i] ) + fabs ( rv1[i] );
973 norm = ( norm > tmp ) ? norm : tmp;
974 }
975
976 // accumulation of right-hand transformations:
977 for ( i = nm1 ; i >= 0 ; --i )
978 {
979 if ( i < nm1 )
980 {
981 if ( g )
982 {
983 for ( j = l ; j < n ; ++j )
984 V[j][i] = ( M[i][j] / M[i][l] ) / g;
985 for ( j = l ; j < n ; ++j ) {
986 for ( s = 0.0 , k = l ; k < n ; ++k )
987 s += M[i][k] * V[k][j];
988 for ( k = l ; k < n ; ++k )
989 V[k][j] += s * V[k][i];
990 }
991 }
992 for ( j = l ; j < n ; ++j )
993 V[i][j] = V[j][i] = 0.0;
994 }
995 V[i][i] = 1.0;
996 g = rv1[i];
997 l = i;
998 }
999
1000 // accumulation of left-hand transformations:
1001 for ( i = ( ( m < n ) ? m : n ) - 1 ; i >= 0 ; --i )
1002 {
1003 l = i + 1;
1004 g = W[i];
1005 for ( j = l ; j < n ; ++j )
1006 M[i][j] = 0.0;
1007 if ( g ) {
1008 g = 1.0 / g;
1009 for ( j = l ; j < n ; ++j )
1010 {
1011 for ( s = 0.0 , k = l ; k < m ; ++k )
1012 s += M[k][i] * M[k][j];
1013 f = ( s / M[i][i] ) * g;
1014 for ( k = i ; k < m ; ++k )
1015 M[k][j] += f * M[k][i];
1016 }
1017 for ( j = i ; j < m ; ++j )
1018 M[j][i] *= g;
1019 }
1020 else
1021 for ( j = i ; j < m ; ++j )
1022 M[j][i] = 0.0;
1023 ++M[i][i];
1024 }
1025
1026 // diagonalization of the bidiagonal form:
1027 for ( k = nm1 ; k >= 0 ; --k )
1028 {
1029 for ( its = 1 ; its <= NITER ; its++ )
1030 {
1031 flag = true;
1032 for ( l = k ; l >= 0 ; l-- )
1033 {
1034 nm = l - 1;
1035 if ( nm < 0 || fabs ( rv1[l]) + norm == norm )
1036 {
1037 flag = false;
1038 break;
1039 }
1040 if ( fabs ( W[nm] ) + norm == norm )
1041 break;
1042 }
1043 if ( flag )
1044 {
1045 c = 0.0;
1046 s = 1.0;
1047 for ( i = l ; i <= k ; i++ )
1048 {
1049 f = s * rv1[i];
1050 rv1[i] = c * rv1[i];
1051 if ( fabs(f) + norm == norm )
1052 break;
1053 g = W[i];
1054
1055 absf = fabs(f);
1056 absg = fabs(g);
1057 h = ( absf > absg ) ?
1058 absf * sqrt ( 1.0 + pow ( absg/absf , 2.0 ) ) :
1059 ( ( absg==0 ) ? 0.0 : absg * sqrt ( 1.0 + pow ( absf/absg , 2.0 ) ) );
1060
1061 W[i] = h;
1062 h = 1.0 / h;
1063 c = g * h;
1064 s = -f * h;
1065 for ( j = 0 ; j < m ; ++j ) {
1066 y = M[j][nm];
1067 z = M[j][ i];
1068 M[j][nm] = y * c + z * s;
1069 M[j][ i] = z * c - y * s;
1070 }
1071 }
1072 }
1073 z = W[k];
1074 if ( l == k) {
1075 if ( z < 0.0 ) {
1076 W[k] = -z;
1077 for ( j = 0 ; j < n ; j++ )
1078 V[j][k] = -V[j][k];
1079 }
1080 break; // this 'break' is always active if k==0
1081 }
1082 if ( its == NITER )
1083 {
1084 error_msg = "SVD_decomposition() error: no convergence in " +
1085 NOMAD::itos ( NITER ) + " iterations";
1086 delete [] rv1;
1087 return false;
1088 }
1089 x = W[l];
1090 nm = k - 1;
1091 y = W[nm];
1092 g = rv1[nm];
1093 h = rv1[k];
1094 f = ( (y-z) * (y+z) + (g-h) * (g+h) ) / ( 2.0 * h * y );
1095
1096 absf = fabs(f);
1097 g = ( absf > 1.0 ) ?
1098 absf * sqrt ( 1.0 + pow ( 1.0/absf , 2.0 ) ) :
1099 sqrt ( 1.0 + pow ( absf , 2.0 ) );
1100
1101 f = ( (x-z) * (x+z) +
1102 h * ( ( y / ( f + ( (f >= 0)? fabs(g) : -fabs(g) ) ) ) - h ) ) / x;
1103 c = s = 1.0;
1104
1105 for ( j = l ; j <= nm ; ++j ) {
1106 i = j + 1;
1107 g = rv1[i];
1108 y = W[i];
1109 h = s * g;
1110 g = c * g;
1111
1112 absf = fabs(f);
1113 absh = fabs(h);
1114 z = ( absf > absh ) ?
1115 absf * sqrt ( 1.0 + pow ( absh/absf , 2.0 ) ) :
1116 ( ( absh==0 ) ? 0.0 : absh * sqrt ( 1.0 + pow ( absf/absh , 2.0 ) ) );
1117
1118 rv1[j] = z;
1119 c = f / z;
1120 s = h / z;
1121 f = x * c + g * s;
1122 g = g * c - x * s;
1123 h = y * s;
1124 y *= c;
1125 for ( jj = 0 ; jj < n ; ++jj )
1126 {
1127 x = V[jj][j];
1128 z = V[jj][i];
1129 V[jj][j] = x * c + z * s;
1130 V[jj][i] = z * c - x * s;
1131 }
1132
1133 absf = fabs(f);
1134 absh = fabs(h);
1135 z = ( absf > absh ) ?
1136 absf * sqrt ( 1.0 + pow ( absh/absf , 2.0 ) ) :
1137 ( ( absh==0 ) ? 0.0 : absh * sqrt ( 1.0 + pow ( absf/absh , 2.0 ) ) );
1138
1139 W[j] = z;
1140
1141 if ( z )
1142 {
1143 z = 1.0 / z;
1144 c = f * z;
1145 s = h * z;
1146 }
1147 f = c * g + s * y;
1148 x = c * y - s * g;
1149 for ( jj = 0 ; jj < m ; ++jj )
1150 {
1151 y = M[jj][j];
1152 z = M[jj][i];
1153 M[jj][j] = y * c + z * s;
1154 M[jj][i] = z * c - y * s;
1155 }
1156 }
1157 rv1[l] = 0.0;
1158 rv1[k] = f;
1159 W [k] = x;
1160 }
1161 }
1162
1163 delete [] rv1;
1164 return true;
1165 }
1166