1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43
44 /** \file
45 \brief Contains definitions of custom data types in ROL.
46 \author Created by D. Ridzal and D. Kouri.
47 */
48
49 #ifndef ROL_TYPES_HPP
50 #define ROL_TYPES_HPP
51
52 #ifdef HAVE_ROL_DEBUG
53 #define ROL_VALIDATE( A ) A
54 #else
55 #define ROL_VALIDATE( A ) /* empty */
56 #endif
57
58 #include <algorithm>
59 #include <complex>
60 #include <exception>
61 #include <string>
62 #include <sstream>
63 #include <limits>
64 #include <type_traits>
65 #include <ROL_stacktrace.hpp>
66 #include "ROL_ScalarTraits.hpp"
67 #include <ROL_Ptr.hpp>
68 #include <ROL_Vector.hpp>
69 #include <ROL_config.h>
70
71 /** \def ROL_NUM_CHECKDERIV_STEPS
72 \brief Number of steps for derivative checks.
73 */
74 #define ROL_NUM_CHECKDERIV_STEPS 13
75
76
77
78 namespace ROL {
79
80 template<class T>
NumberToString(T Number)81 std::string NumberToString( T Number )
82 {
83 std::ostringstream ss;
84 ss << Number;
85 return ss.str();
86 }
87
88 /** \brief Platform-dependent machine epsilon.
89 */
90 template<class Real>
ROL_EPSILON(void)91 inline Real ROL_EPSILON(void) { return std::abs(ROL::ScalarTraits<Real>::eps()); }
92 //static const Real ROL_EPSILON<Real>() = std::abs(ROL::ScalarTraits<Real>::eps());
93
94 /** \brief Tolerance for various equality tests.
95 */
96 template<class Real>
ROL_THRESHOLD(void)97 inline Real ROL_THRESHOLD(void) { return 10.0 * ROL_EPSILON<Real>(); }
98
99 /** \brief Platform-dependent maximum double.
100 */
101 template<class Real>
ROL_OVERFLOW(void)102 inline Real ROL_OVERFLOW(void) { return std::abs(ROL::ScalarTraits<Real>::rmax()); }
103
104 template<class Real>
ROL_INF(void)105 inline Real ROL_INF(void) { return 0.1*ROL_OVERFLOW<Real>(); }
106
107 template<class Real>
ROL_NINF(void)108 inline Real ROL_NINF(void) { return -ROL_INF<Real>(); }
109
110 /** \brief Platform-dependent minimum double.
111 */
112 template<class Real>
ROL_UNDERFLOW(void)113 inline Real ROL_UNDERFLOW(void) { return std::abs(ROL::ScalarTraits<Real>::rmin()); }
114
115 /** \brief Enum for algorithm termination.
116 */
117 enum EExitStatus {
118 EXITSTATUS_CONVERGED = 0,
119 EXITSTATUS_MAXITER,
120 EXITSTATUS_STEPTOL,
121 EXITSTATUS_NAN,
122 EXITSTATUS_USERDEFINED,
123 EXITSTATUS_LAST
124 };
125
EExitStatusToString(EExitStatus tr)126 inline std::string EExitStatusToString(EExitStatus tr) {
127 std::string retString;
128 switch(tr) {
129 case EXITSTATUS_CONVERGED: retString = "Converged"; break;
130 case EXITSTATUS_MAXITER: retString = "Iteration Limit Exceeded"; break;
131 case EXITSTATUS_STEPTOL: retString = "Step Tolerance Met"; break;
132 case EXITSTATUS_NAN: retString = "Step and/or Gradient Returned NaN"; break;
133 case EXITSTATUS_USERDEFINED: retString = "User Defined"; break;
134 case EXITSTATUS_LAST: retString = "Last Type (Dummy)"; break;
135 default: retString = "INVALID EExitStatus";
136 }
137 return retString;
138 }
139
140 /** \brief State for algorithm class. Will be used for restarts.
141 */
142 template<class Real>
143 struct AlgorithmState {
144 int iter;
145 int minIter;
146 int nfval;
147 int ncval;
148 int ngrad;
149 Real value;
150 Real minValue;
151 Real gnorm;
152 Real cnorm;
153 Real snorm;
154 Real aggregateGradientNorm;
155 Real aggregateModelError;
156 bool flag;
157 ROL::Ptr<Vector<Real> > iterateVec;
158 ROL::Ptr<Vector<Real> > lagmultVec;
159 ROL::Ptr<Vector<Real> > minIterVec;
160 EExitStatus statusFlag;
161
AlgorithmStateROL::AlgorithmState162 AlgorithmState(void) : iter(0), minIter(0), nfval(0), ngrad(0), value(0), minValue(0),
163 gnorm(std::numeric_limits<Real>::max()),
164 cnorm(std::numeric_limits<Real>::max()),
165 snorm(std::numeric_limits<Real>::max()),
166 aggregateGradientNorm(std::numeric_limits<Real>::max()),
167 aggregateModelError(std::numeric_limits<Real>::max()),
168 flag(false),
169 iterateVec(ROL::nullPtr), lagmultVec(ROL::nullPtr), minIterVec(ROL::nullPtr),
170 statusFlag(EXITSTATUS_LAST) {}
171
resetROL::AlgorithmState172 void reset(void) {
173 iter = 0;
174 minIter = 0;
175 nfval = 0;
176 ncval = 0;
177 ngrad = 0;
178 value = ROL_INF<Real>();
179 minValue = ROL_INF<Real>();
180 gnorm = ROL_INF<Real>();
181 cnorm = ROL_INF<Real>();
182 snorm = ROL_INF<Real>();
183 aggregateGradientNorm = ROL_INF<Real>();
184 aggregateModelError = ROL_INF<Real>();
185 flag = false;
186 if (iterateVec != ROL::nullPtr) {
187 iterateVec->zero();
188 }
189 if (lagmultVec != ROL::nullPtr) {
190 lagmultVec->zero();
191 }
192 if (minIterVec != ROL::nullPtr) {
193 minIterVec->zero();
194 }
195 }
196 };
197
198 /** \brief State for step class. Will be used for restarts.
199 */
200 template<class Real>
201 struct StepState {
202 ROL::Ptr<Vector<Real> > gradientVec;
203 ROL::Ptr<Vector<Real> > descentVec;
204 ROL::Ptr<Vector<Real> > constraintVec;
205 int nfval;
206 int ngrad;
207 Real searchSize; // line search parameter (alpha) or trust-region radius (delta)
208 int flag; // Was step successful?
209 int SPiter; // Subproblem iteration count
210 int SPflag; // Subproblem termination flag
211
StepStateROL::StepState212 StepState(void) : gradientVec(ROL::nullPtr),
213 descentVec(ROL::nullPtr),
214 constraintVec(ROL::nullPtr),
215 nfval(0),
216 ngrad(0),
217 searchSize(0),
218 flag(0),
219 SPiter(0),
220 SPflag(0) {}
221
resetROL::StepState222 void reset(const Real searchSizeInput = 1.0) {
223 if (gradientVec != ROL::nullPtr) {
224 gradientVec->zero();
225 }
226 if (descentVec != ROL::nullPtr) {
227 descentVec->zero();
228 }
229 if (constraintVec != ROL::nullPtr) {
230 constraintVec->zero();
231 }
232 nfval = 0;
233 ngrad = 0;
234 searchSize = searchSizeInput;
235 flag = 0;
236 SPiter = 0;
237 SPflag = 0;
238 }
239 };
240
241 struct removeSpecialCharacters {
operator ()ROL::removeSpecialCharacters242 bool operator()(char c) {
243 return (c ==' ' || c =='-' || c == '(' || c == ')' || c=='\'' || c=='\r' || c=='\n' || c=='\t');
244 }
245 };
246
removeStringFormat(std::string s)247 inline std::string removeStringFormat( std::string s ) {
248 std::string output = s;
249 output.erase( std::remove_if( output.begin(), output.end(), removeSpecialCharacters()), output.end() );
250 std::transform( output.begin(), output.end(), output.begin(), ::tolower );
251 return output;
252 }
253
254 // Types of optimization problem
255 enum EProblem {
256 TYPE_U = 0,
257 TYPE_B,
258 TYPE_E,
259 TYPE_EB,
260 TYPE_LAST
261 };
262
263 /** \enum ROL::EStep
264 \brief Enumeration of step types.
265
266 \arg AUGMENTEDLAGRANGIAN describe
267 \arg BUNDLE describe
268 \arg COMPOSITESTEP describe
269 \arg LINESEARCH describe
270 \arg MOREAUYOSIDAPENALTY describe
271 \arg PRIMALDUALACTIVESET describe
272 \arg TRUSTREGION describe
273 */
274 enum EStep{
275 STEP_AUGMENTEDLAGRANGIAN = 0,
276 STEP_BUNDLE,
277 STEP_COMPOSITESTEP,
278 STEP_LINESEARCH,
279 STEP_MOREAUYOSIDAPENALTY,
280 STEP_PRIMALDUALACTIVESET,
281 STEP_TRUSTREGION,
282 STEP_INTERIORPOINT,
283 STEP_FLETCHER,
284 STEP_LAST
285 };
286
EStepToString(EStep tr)287 inline std::string EStepToString(EStep tr) {
288 std::string retString;
289 switch(tr) {
290 case STEP_AUGMENTEDLAGRANGIAN: retString = "Augmented Lagrangian"; break;
291 case STEP_BUNDLE: retString = "Bundle"; break;
292 case STEP_COMPOSITESTEP: retString = "Composite Step"; break;
293 case STEP_LINESEARCH: retString = "Line Search"; break;
294 case STEP_MOREAUYOSIDAPENALTY: retString = "Moreau-Yosida Penalty"; break;
295 case STEP_PRIMALDUALACTIVESET: retString = "Primal Dual Active Set"; break;
296 case STEP_TRUSTREGION: retString = "Trust Region"; break;
297 case STEP_INTERIORPOINT: retString = "Interior Point"; break;
298 case STEP_FLETCHER: retString = "Fletcher"; break;
299 case STEP_LAST: retString = "Last Type (Dummy)"; break;
300 default: retString = "INVALID EStep";
301 }
302 return retString;
303 }
304
isCompatibleStep(EProblem p,EStep s)305 inline bool isCompatibleStep( EProblem p, EStep s ) {
306 bool comp = false;
307 switch(p) {
308
309 case TYPE_U: comp = ( (s == STEP_LINESEARCH) ||
310 (s == STEP_TRUSTREGION) ||
311 (s == STEP_BUNDLE) );
312 break;
313
314 case TYPE_B: comp = ( (s == STEP_LINESEARCH) ||
315 (s == STEP_TRUSTREGION) ||
316 (s == STEP_MOREAUYOSIDAPENALTY) ||
317 (s == STEP_PRIMALDUALACTIVESET) ||
318 (s == STEP_INTERIORPOINT) );
319 break;
320
321 case TYPE_E: comp = ( (s == STEP_COMPOSITESTEP) ||
322 (s == STEP_AUGMENTEDLAGRANGIAN) ||
323 (s == STEP_FLETCHER) );
324 break;
325
326 case TYPE_EB: comp = ( (s == STEP_AUGMENTEDLAGRANGIAN) ||
327 (s == STEP_MOREAUYOSIDAPENALTY) ||
328 (s == STEP_INTERIORPOINT) ||
329 (s == STEP_FLETCHER) );
330 break;
331
332 case TYPE_LAST: comp = false; break;
333 default: comp = false;
334 }
335 return comp;
336 }
337
EProblemToString(EProblem p)338 inline std::string EProblemToString( EProblem p ) {
339 std::string retString;
340 switch(p) {
341 case TYPE_U: retString = "Type-U"; break;
342 case TYPE_E: retString = "Type-E"; break;
343 case TYPE_B: retString = "Type-B"; break;
344 case TYPE_EB: retString = "Type-EB"; break;
345 case TYPE_LAST: retString = "Type-Last (Dummy)"; break;
346 default: retString = "Invalid EProblem";
347 }
348 return retString;
349 }
350
351
352 /** \brief Verifies validity of a TrustRegion enum.
353
354 \param tr [in] - enum of the TrustRegion
355 \return 1 if the argument is a valid TrustRegion; 0 otherwise.
356 */
isValidStep(EStep ls)357 inline int isValidStep(EStep ls) {
358 return( (ls == STEP_AUGMENTEDLAGRANGIAN) ||
359 (ls == STEP_BUNDLE) ||
360 (ls == STEP_COMPOSITESTEP) ||
361 (ls == STEP_LINESEARCH) ||
362 (ls == STEP_MOREAUYOSIDAPENALTY) ||
363 (ls == STEP_PRIMALDUALACTIVESET) ||
364 (ls == STEP_TRUSTREGION) ||
365 (ls == STEP_INTERIORPOINT) ||
366 (ls == STEP_FLETCHER) ) ;
367 }
368
operator ++(EStep & type)369 inline EStep & operator++(EStep &type) {
370 return type = static_cast<EStep>(type+1);
371 }
372
operator ++(EStep & type,int)373 inline EStep operator++(EStep &type, int) {
374 EStep oldval = type;
375 ++type;
376 return oldval;
377 }
378
operator --(EStep & type)379 inline EStep & operator--(EStep &type) {
380 return type = static_cast<EStep>(type-1);
381 }
382
operator --(EStep & type,int)383 inline EStep operator--(EStep &type, int) {
384 EStep oldval = type;
385 --type;
386 return oldval;
387 }
388
StringToEStep(std::string s)389 inline EStep StringToEStep(std::string s) {
390 s = removeStringFormat(s);
391 for ( EStep st = STEP_AUGMENTEDLAGRANGIAN; st < STEP_LAST; ++st ) {
392 if ( !s.compare(removeStringFormat(EStepToString(st))) ) {
393 return st;
394 }
395 }
396 return STEP_LAST;
397 }
398
399 /** \enum ROL::EDescent
400 \brief Enumeration of descent direction types.
401
402 \arg STEEPEST describe
403 \arg NONLINEARCG describe
404 \arg SECANT describe
405 \arg NEWTON describe
406 \arg NEWTONKRYLOV describe
407 \arg SECANTPRECOND describe
408 */
409 enum EDescent{
410 DESCENT_STEEPEST = 0,
411 DESCENT_NONLINEARCG,
412 DESCENT_SECANT,
413 DESCENT_NEWTON,
414 DESCENT_NEWTONKRYLOV,
415 DESCENT_LAST
416 };
417
EDescentToString(EDescent tr)418 inline std::string EDescentToString(EDescent tr) {
419 std::string retString;
420 switch(tr) {
421 case DESCENT_STEEPEST: retString = "Steepest Descent"; break;
422 case DESCENT_NONLINEARCG: retString = "Nonlinear CG"; break;
423 case DESCENT_SECANT: retString = "Quasi-Newton Method"; break;
424 case DESCENT_NEWTON: retString = "Newton's Method"; break;
425 case DESCENT_NEWTONKRYLOV: retString = "Newton-Krylov"; break;
426 case DESCENT_LAST: retString = "Last Type (Dummy)"; break;
427 default: retString = "INVALID ESecant";
428 }
429 return retString;
430 }
431
432 /** \brief Verifies validity of a Secant enum.
433
434 \param tr [in] - enum of the Secant
435 \return 1 if the argument is a valid Secant; 0 otherwise.
436 */
isValidDescent(EDescent d)437 inline int isValidDescent(EDescent d){
438 return( (d == DESCENT_STEEPEST) ||
439 (d == DESCENT_NONLINEARCG) ||
440 (d == DESCENT_SECANT) ||
441 (d == DESCENT_NEWTON) ||
442 (d == DESCENT_NEWTONKRYLOV)
443 );
444 }
445
operator ++(EDescent & type)446 inline EDescent & operator++(EDescent &type) {
447 return type = static_cast<EDescent>(type+1);
448 }
449
operator ++(EDescent & type,int)450 inline EDescent operator++(EDescent &type, int) {
451 EDescent oldval = type;
452 ++type;
453 return oldval;
454 }
455
operator --(EDescent & type)456 inline EDescent & operator--(EDescent &type) {
457 return type = static_cast<EDescent>(type-1);
458 }
459
operator --(EDescent & type,int)460 inline EDescent operator--(EDescent &type, int) {
461 EDescent oldval = type;
462 --type;
463 return oldval;
464 }
465
StringToEDescent(std::string s)466 inline EDescent StringToEDescent(std::string s) {
467 s = removeStringFormat(s);
468 for ( EDescent des = DESCENT_STEEPEST; des < DESCENT_LAST; des++ ) {
469 if ( !s.compare(removeStringFormat(EDescentToString(des))) ) {
470 return des;
471 }
472 }
473 return DESCENT_SECANT;
474 }
475
476 /** \enum ROL::ESecant
477 \brief Enumeration of secant update algorithms.
478
479 \arg LBFGS describe
480 \arg LDFP describe
481 \arg LSR1 describe
482 \arg BARZILAIBORWEIN describe
483 */
484 enum ESecant{
485 SECANT_LBFGS = 0,
486 SECANT_LDFP,
487 SECANT_LSR1,
488 SECANT_BARZILAIBORWEIN,
489 SECANT_USERDEFINED,
490 SECANT_LAST
491 };
492
ESecantToString(ESecant tr)493 inline std::string ESecantToString(ESecant tr) {
494 std::string retString;
495 switch(tr) {
496 case SECANT_LBFGS: retString = "Limited-Memory BFGS"; break;
497 case SECANT_LDFP: retString = "Limited-Memory DFP"; break;
498 case SECANT_LSR1: retString = "Limited-Memory SR1"; break;
499 case SECANT_BARZILAIBORWEIN: retString = "Barzilai-Borwein"; break;
500 case SECANT_USERDEFINED: retString = "User-Defined"; break;
501 case SECANT_LAST: retString = "Last Type (Dummy)"; break;
502 default: retString = "INVALID ESecant";
503 }
504 return retString;
505 }
506
507 /** \brief Verifies validity of a Secant enum.
508
509 \param tr [in] - enum of the Secant
510 \return 1 if the argument is a valid Secant; 0 otherwise.
511 */
isValidSecant(ESecant s)512 inline int isValidSecant(ESecant s) {
513 return( (s == SECANT_LBFGS) ||
514 (s == SECANT_LDFP) ||
515 (s == SECANT_LSR1) ||
516 (s == SECANT_BARZILAIBORWEIN) ||
517 (s == SECANT_USERDEFINED)
518 );
519 }
520
operator ++(ESecant & type)521 inline ESecant & operator++(ESecant &type) {
522 return type = static_cast<ESecant>(type+1);
523 }
524
operator ++(ESecant & type,int)525 inline ESecant operator++(ESecant &type, int) {
526 ESecant oldval = type;
527 ++type;
528 return oldval;
529 }
530
operator --(ESecant & type)531 inline ESecant & operator--(ESecant &type) {
532 return type = static_cast<ESecant>(type-1);
533 }
534
operator --(ESecant & type,int)535 inline ESecant operator--(ESecant &type, int) {
536 ESecant oldval = type;
537 --type;
538 return oldval;
539 }
540
StringToESecant(std::string s)541 inline ESecant StringToESecant(std::string s) {
542 s = removeStringFormat(s);
543 for ( ESecant sec = SECANT_LBFGS; sec < SECANT_LAST; sec++ ) {
544 if ( !s.compare(removeStringFormat(ESecantToString(sec))) ) {
545 return sec;
546 }
547 }
548 return SECANT_LBFGS;
549 }
550
551 /** \enum ROL::ENonlinearCG
552 \brief Enumeration of nonlinear CG algorithms.
553
554 \arg HESTENES_STIEFEL \f$ \frac{g_{k+1}^\top y_k}{d_k^\top y_k } \f$
555 \arg FLETCHER_REEVES \f$ \frac{\|g_{k+1}\|^2}{\|g_k\|^2} \f$
556 \arg DANIEL \f$ \frac{g_{k+1}^\top \nabla^2 f(x_k) d_k}{d_k^\top \nabla^2 f(x_k) d_k} \f$
557 \arg POLAK_RIBIERE \f$ \frac{g_{k+1}^\top y_k}{\|g_k\|^2} \f$
558 \arg FLETCHER_CONJDESC \f$ -\frac{\|g_{k+1}\|^2}{d_k^\top g_k} \f$
559 \arg LIU_STOREY \f$ -\frac{g_k^\top y_{k-1} }{d_{k-1}^\top g_{k-1} \f$
560 \arg DAI_YUAN \f$ \frac{\|g_{k+1}\|^2}{d_k^\top y_k} \f$
561 \arg HAGER_ZHANG \f$ \frac{g_{k+1}^\top y_k}{d_k^\top y_k} - 2 \frac{\|y_k\|^2}{d_k^\top y_k} \frac{g_{k+1}^\top d_k}{d_k^\top y_k} \f$
562 \arg OREN_LUENBERGER \f$ \frac{g_{k+1}^\top y_k}{d_k^\top y_k} - \frac{\|y_k\|^2}{d_k^\top y_k} \frac{g_{k+1}^\top d_k}{d_k^\top y_k} \f$
563 */
564 enum ENonlinearCG{
565 NONLINEARCG_HESTENES_STIEFEL = 0,
566 NONLINEARCG_FLETCHER_REEVES,
567 NONLINEARCG_DANIEL,
568 NONLINEARCG_POLAK_RIBIERE,
569 NONLINEARCG_FLETCHER_CONJDESC,
570 NONLINEARCG_LIU_STOREY,
571 NONLINEARCG_DAI_YUAN,
572 NONLINEARCG_HAGER_ZHANG,
573 NONLINEARCG_OREN_LUENBERGER,
574 NONLINEARCG_USERDEFINED,
575 NONLINEARCG_LAST
576 };
577
ENonlinearCGToString(ENonlinearCG tr)578 inline std::string ENonlinearCGToString(ENonlinearCG tr) {
579 std::string retString;
580 switch(tr) {
581 case NONLINEARCG_HESTENES_STIEFEL: retString = "Hestenes-Stiefel"; break;
582 case NONLINEARCG_FLETCHER_REEVES: retString = "Fletcher-Reeves"; break;
583 case NONLINEARCG_DANIEL: retString = "Daniel (uses Hessian)"; break;
584 case NONLINEARCG_POLAK_RIBIERE: retString = "Polak-Ribiere"; break;
585 case NONLINEARCG_FLETCHER_CONJDESC: retString = "Fletcher Conjugate Descent"; break;
586 case NONLINEARCG_LIU_STOREY: retString = "Liu-Storey"; break;
587 case NONLINEARCG_DAI_YUAN: retString = "Dai-Yuan"; break;
588 case NONLINEARCG_HAGER_ZHANG: retString = "Hager-Zhang"; break;
589 case NONLINEARCG_OREN_LUENBERGER: retString = "Oren-Luenberger"; break;
590 case NONLINEARCG_USERDEFINED: retString = "User Defined"; break;
591 case NONLINEARCG_LAST: retString = "Last Type (Dummy)"; break;
592 default: retString = "INVALID ENonlinearCG";
593 }
594 return retString;
595 }
596
597 /** \brief Verifies validity of a NonlinearCG enum.
598
599 \param tr [in] - enum of the NonlinearCG
600 \return 1 if the argument is a valid NonlinearCG; 0 otherwise.
601 */
isValidNonlinearCG(ENonlinearCG s)602 inline int isValidNonlinearCG(ENonlinearCG s) {
603 return( (s == NONLINEARCG_HESTENES_STIEFEL) ||
604 (s == NONLINEARCG_FLETCHER_REEVES) ||
605 (s == NONLINEARCG_DANIEL) ||
606 (s == NONLINEARCG_POLAK_RIBIERE) ||
607 (s == NONLINEARCG_FLETCHER_CONJDESC) ||
608 (s == NONLINEARCG_LIU_STOREY) ||
609 (s == NONLINEARCG_DAI_YUAN) ||
610 (s == NONLINEARCG_HAGER_ZHANG) ||
611 (s == NONLINEARCG_OREN_LUENBERGER) ||
612 (s == NONLINEARCG_USERDEFINED)
613 );
614 }
615
operator ++(ENonlinearCG & type)616 inline ENonlinearCG & operator++(ENonlinearCG &type) {
617 return type = static_cast<ENonlinearCG>(type+1);
618 }
619
operator ++(ENonlinearCG & type,int)620 inline ENonlinearCG operator++(ENonlinearCG &type, int) {
621 ENonlinearCG oldval = type;
622 ++type;
623 return oldval;
624 }
625
operator --(ENonlinearCG & type)626 inline ENonlinearCG & operator--(ENonlinearCG &type) {
627 return type = static_cast<ENonlinearCG>(type-1);
628 }
629
operator --(ENonlinearCG & type,int)630 inline ENonlinearCG operator--(ENonlinearCG &type, int) {
631 ENonlinearCG oldval = type;
632 --type;
633 return oldval;
634 }
635
StringToENonlinearCG(std::string s)636 inline ENonlinearCG StringToENonlinearCG(std::string s) {
637 s = removeStringFormat(s);
638 for ( ENonlinearCG nlcg = NONLINEARCG_HESTENES_STIEFEL; nlcg < NONLINEARCG_LAST; nlcg++ ) {
639 if ( !s.compare(removeStringFormat(ENonlinearCGToString(nlcg))) ) {
640 return nlcg;
641 }
642 }
643 return NONLINEARCG_HESTENES_STIEFEL;
644 }
645
646 /** \enum ROL::ELineSearch
647 \brief Enumeration of line-search types.
648
649 \arg BACKTRACKING describe
650 \arg BISECTION describe
651 \arg GOLDENSECTION describe
652 \arg CUBICINTERP describe
653 \arg BRENTS describe
654 \arg USERDEFINED describe
655 */
656 enum ELineSearch{
657 LINESEARCH_ITERATIONSCALING = 0,
658 LINESEARCH_PATHBASEDTARGETLEVEL,
659 LINESEARCH_BACKTRACKING,
660 LINESEARCH_BISECTION,
661 LINESEARCH_GOLDENSECTION,
662 LINESEARCH_CUBICINTERP,
663 LINESEARCH_BRENTS,
664 LINESEARCH_USERDEFINED,
665 LINESEARCH_LAST
666 };
667
ELineSearchToString(ELineSearch ls)668 inline std::string ELineSearchToString(ELineSearch ls) {
669 std::string retString;
670 switch(ls) {
671 case LINESEARCH_ITERATIONSCALING: retString = "Iteration Scaling"; break;
672 case LINESEARCH_PATHBASEDTARGETLEVEL: retString = "Path-Based Target Level"; break;
673 case LINESEARCH_BACKTRACKING: retString = "Backtracking"; break;
674 case LINESEARCH_BISECTION: retString = "Bisection"; break;
675 case LINESEARCH_GOLDENSECTION: retString = "Golden Section"; break;
676 case LINESEARCH_CUBICINTERP: retString = "Cubic Interpolation"; break;
677 case LINESEARCH_BRENTS: retString = "Brent's"; break;
678 case LINESEARCH_USERDEFINED: retString = "User Defined"; break;
679 case LINESEARCH_LAST: retString = "Last Type (Dummy)"; break;
680 default: retString = "INVALID ELineSearch";
681 }
682 return retString;
683 }
684
685 /** \brief Verifies validity of a LineSearch enum.
686
687 \param ls [in] - enum of the linesearch
688 \return 1 if the argument is a valid linesearch; 0 otherwise.
689 */
isValidLineSearch(ELineSearch ls)690 inline int isValidLineSearch(ELineSearch ls){
691 return( (ls == LINESEARCH_BACKTRACKING) ||
692 (ls == LINESEARCH_ITERATIONSCALING) ||
693 (ls == LINESEARCH_PATHBASEDTARGETLEVEL) ||
694 (ls == LINESEARCH_BISECTION) ||
695 (ls == LINESEARCH_GOLDENSECTION) ||
696 (ls == LINESEARCH_CUBICINTERP) ||
697 (ls == LINESEARCH_BRENTS) ||
698 (ls == LINESEARCH_USERDEFINED)
699 );
700 }
701
operator ++(ELineSearch & type)702 inline ELineSearch & operator++(ELineSearch &type) {
703 return type = static_cast<ELineSearch>(type+1);
704 }
705
operator ++(ELineSearch & type,int)706 inline ELineSearch operator++(ELineSearch &type, int) {
707 ELineSearch oldval = type;
708 ++type;
709 return oldval;
710 }
711
operator --(ELineSearch & type)712 inline ELineSearch & operator--(ELineSearch &type) {
713 return type = static_cast<ELineSearch>(type-1);
714 }
715
operator --(ELineSearch & type,int)716 inline ELineSearch operator--(ELineSearch &type, int) {
717 ELineSearch oldval = type;
718 --type;
719 return oldval;
720 }
721
StringToELineSearch(std::string s)722 inline ELineSearch StringToELineSearch(std::string s) {
723 s = removeStringFormat(s);
724 for ( ELineSearch ls = LINESEARCH_ITERATIONSCALING; ls < LINESEARCH_LAST; ls++ ) {
725 if ( !s.compare(removeStringFormat(ELineSearchToString(ls))) ) {
726 return ls;
727 }
728 }
729 return LINESEARCH_ITERATIONSCALING;
730 }
731
732 /** \enum ROL::ECurvatureCondition
733 \brief Enumeration of line-search curvature conditions.
734
735 \arg WOLFE describe
736 \arg STRONGWOLFE describe
737 \arg GOLDSTEIN describe
738 */
739 enum ECurvatureCondition{
740 CURVATURECONDITION_WOLFE = 0,
741 CURVATURECONDITION_STRONGWOLFE,
742 CURVATURECONDITION_GENERALIZEDWOLFE,
743 CURVATURECONDITION_APPROXIMATEWOLFE,
744 CURVATURECONDITION_GOLDSTEIN,
745 CURVATURECONDITION_NULL,
746 CURVATURECONDITION_LAST
747 };
748
ECurvatureConditionToString(ECurvatureCondition ls)749 inline std::string ECurvatureConditionToString(ECurvatureCondition ls) {
750 std::string retString;
751 switch(ls) {
752 case CURVATURECONDITION_WOLFE: retString = "Wolfe Conditions"; break;
753 case CURVATURECONDITION_STRONGWOLFE: retString = "Strong Wolfe Conditions"; break;
754 case CURVATURECONDITION_GENERALIZEDWOLFE: retString = "Generalized Wolfe Conditions"; break;
755 case CURVATURECONDITION_APPROXIMATEWOLFE: retString = "Approximate Wolfe Conditions"; break;
756 case CURVATURECONDITION_GOLDSTEIN: retString = "Goldstein Conditions"; break;
757 case CURVATURECONDITION_NULL: retString = "Null Curvature Condition"; break;
758 case CURVATURECONDITION_LAST: retString = "Last Type (Dummy)"; break;
759 default: retString = "INVALID ECurvatureCondition";
760 }
761 return retString;
762 }
763
764 /** \brief Verifies validity of a CurvatureCondition enum.
765
766 \param ls [in] - enum of the Curvature Conditions
767 \return 1 if the argument is a valid curvature condition; 0 otherwise.
768 */
isValidCurvatureCondition(ECurvatureCondition ls)769 inline int isValidCurvatureCondition(ECurvatureCondition ls){
770 return( (ls == CURVATURECONDITION_WOLFE) ||
771 (ls == CURVATURECONDITION_STRONGWOLFE) ||
772 (ls == CURVATURECONDITION_GENERALIZEDWOLFE) ||
773 (ls == CURVATURECONDITION_APPROXIMATEWOLFE) ||
774 (ls == CURVATURECONDITION_GOLDSTEIN) ||
775 (ls == CURVATURECONDITION_NULL)
776 );
777 }
778
operator ++(ECurvatureCondition & type)779 inline ECurvatureCondition & operator++(ECurvatureCondition &type) {
780 return type = static_cast<ECurvatureCondition>(type+1);
781 }
782
operator ++(ECurvatureCondition & type,int)783 inline ECurvatureCondition operator++(ECurvatureCondition &type, int) {
784 ECurvatureCondition oldval = type;
785 ++type;
786 return oldval;
787 }
788
operator --(ECurvatureCondition & type)789 inline ECurvatureCondition & operator--(ECurvatureCondition &type) {
790 return type = static_cast<ECurvatureCondition>(type-1);
791 }
792
operator --(ECurvatureCondition & type,int)793 inline ECurvatureCondition operator--(ECurvatureCondition &type, int) {
794 ECurvatureCondition oldval = type;
795 --type;
796 return oldval;
797 }
798
StringToECurvatureCondition(std::string s)799 inline ECurvatureCondition StringToECurvatureCondition(std::string s) {
800 s = removeStringFormat(s);
801 for ( ECurvatureCondition cc = CURVATURECONDITION_WOLFE; cc < CURVATURECONDITION_LAST; cc++ ) {
802 if ( !s.compare(removeStringFormat(ECurvatureConditionToString(cc))) ) {
803 return cc;
804 }
805 }
806 return CURVATURECONDITION_WOLFE;
807 }
808
809 /** \enum ROL::ECGFlag
810 \brief Enumation of flags used by conjugate gradient methods.
811
812 \arg CG_FLAG_SUCCESS Residual Tolerance Met
813 \arg CG_FLAG_ITEREXCEED Iteration Limit Exceeded
814 \arg CG_FLAG_NEGCURVE Negative Curvature Detected
815 \arh CG_FLAG_TRRADEX Trust-Region Radius Exceeded
816 \arh CG_FLAG_ZERORHS Initiali Right Hand Side is Zero
817
818 */
819 enum ECGFlag {
820 CG_FLAG_SUCCESS = 0,
821 CG_FLAG_ITEREXCEED,
822 CG_FLAG_NEGCURVE,
823 CG_FLAG_TRRADEX,
824 CG_FLAG_ZERORHS,
825 CG_FLAG_UNDEFINED
826 };
827
828
ECGFlagToString(ECGFlag cgf)829 inline std::string ECGFlagToString(ECGFlag cgf) {
830 std::string retString;
831 switch(cgf) {
832 case CG_FLAG_SUCCESS:
833 retString = "Residual tolerance met";
834 break;
835 case CG_FLAG_ITEREXCEED:
836 retString = "Iteration limit exceeded";
837 break;
838 case CG_FLAG_NEGCURVE:
839 retString = "Negative curvature detected";
840 break;
841 case CG_FLAG_TRRADEX:
842 retString = "Trust-Region radius exceeded";
843 break;
844 case CG_FLAG_ZERORHS:
845 retString = "Initial right hand side is zero";
846 break;
847 default:
848 retString = "INVALID ECGFlag";
849 }
850 return retString;
851 }
852
853
854
855 // For use in gradient and Hessian checks
856 namespace Finite_Difference_Arrays {
857
858 // Finite difference steps in axpy form
859 const int shifts[4][4] = { { 1, 0, 0, 0 }, // First order
860 { -1, 2, 0, 0 }, // Second order
861 { -1, 2, 1, 0 }, // Third order
862 { -1, -1, 3, 1 } // Fourth order
863 };
864
865 // Finite difference weights
866 const double weights[4][5] = { { -1.0, 1.0, 0.0, 0.0, 0.0 }, // First order
867 { 0.0, -1.0/2.0, 1.0/2.0, 0.0, 0.0 }, // Second order
868 { -1.0/2.0, -1.0/3.0, 1.0, -1.0/6.0, 0.0 }, // Third order
869 { 0.0, -2.0/3.0, 1.0/12.0, 2.0/3.0, -1.0/12.0 } // Fourth order
870 };
871
872 }
873
874
875 // Generic conversion from Element type to Real type
876 template<class Real, class Element>
877 struct TypeCaster {
ElementToRealROL::TypeCaster878 static Real ElementToReal( const Element &val ) {
879 return Real(0);
880 }
881 };
882
883 // Partially specialize for complex<Real>
884 template<class Real>
885 struct TypeCaster<Real, std::complex<Real> > {
ElementToRealROL::TypeCaster886 static Real ElementToReal( const std::complex<Real> &val ) {
887 return val.real();
888 }
889 };
890
891 // Fully specialize for double,float
892 template<>
893 struct TypeCaster<double,float> {
ElementToRealROL::TypeCaster894 static double ElementToReal( const float &val ) {
895 return static_cast<double>(val);
896 }
897 };
898
899 // Cast from Element type to Real type
900 template<class Element, class Real>
rol_cast(const Element & val)901 Real rol_cast(const Element &val) {
902 return TypeCaster<Real,Element>::ElementToReal(val);
903 }
904
905
906
907
908
909
910 namespace Exception {
911
912 class NotImplemented : public std::logic_error {
913 public:
NotImplemented(const std::string & what_arg)914 NotImplemented( const std::string& what_arg ) :
915 std::logic_error(what_arg) {}
916
917
918 }; // class NotImplemented
919
920
921 #if __cplusplus >= 201402L // using C++14
922
923 using std::enable_if_t;
924
925 #else // No C++14
926
927 template<bool B, class T=void>
928 using enable_if_t = typename std::enable_if<B,T>::type;
929
930 #endif
931
932
933
934
935
936 } // namespace Exception
937
938
939 } // namespace ROL
940
941
942 /*! \mainpage %ROL Documentation (Development Version)
943
944 \image html rol.png "Rapid Optimization Library" width=1in
945 \image latex rol.pdf "Rapid Optimization Library" width=1in
946
947 \section intro_sec Introduction
948
949 Rapid Optimization Library (%ROL) is a C++ package for large-scale
950 optimization.
951 It is used for the solution of optimal design, optimal control and
952 inverse problems in large-scale engineering applications.
953 Other uses include mesh optimization and image processing.
954
955
956 \section overview_sec Overview
957
958 %ROL aims to combine flexibility, efficiency and robustness. Key features:
959
960 \li Matrix-free application programming interfaces (APIs) ---enable direct
961 use of application data structures and memory spaces, linear solvers,
962 nonlinear solvers and preconditioners.
963 \li State-of-the-art algorithms for unconstrained optimization,
964 constrained optimization and optimization under uncertainty ---enable
965 inexact and adaptive function evaluations and iterative linear
966 system solves.
967 \li Special APIs for simulation-based optimization ---enable a
968 streamlined embedding into engineering applications, rigorous
969 implementation verification and efficient use.
970 \li Modular interfaces throughout the optimization process ---enable custom
971 and user-defined algorithms, stopping criteria, hierarchies of
972 algorithms, and selective use of a variety of tools and components.
973
974 For a detailed description of user interfaces and algorithms, see the
975 presentations ROL-Trilinos-xx.x.pptx (or .pdf) in the doc/presentations
976 directory.
977
978 To start using %ROL, including all its advanced algorithms and features,
979 jump to the <a href="modules.html">Modules</a> page.
980
981 For a basic example, see below.
982
983 \section quickstart_sec Quick Start
984
985 The Rosenbrock example (rol/example/rosenbrock/example_01.cpp) demonstrates
986 the basic use of %ROL.
987 It amounts to six steps:
988
989 \subsection vector_qs_sec Step 1: Implement linear algebra / vector interface.
990 --- or try one of the provided implementations, such as ROL::StdVector in rol/vector.
991
992 ~~~{.hpp}
993 ROL::Vector
994 ~~~
995
996 \subsection objective_qs_sec Step 2: Implement objective function interface.
997 --- or try one of the provided functions, such as @b ROL::ZOO::Objective_Rosenbrock in rol/zoo.
998
999 \code
1000 ROL::Objective
1001 \endcode
1002
1003 \subsection step_qs_sec Step 3: Choose optimization algorithm.
1004 --- with @b ROL::ParameterList settings in the variable @b parlist.
1005
1006 \code
1007 ROL::Algorithm<RealT> algo("Line Search",parlist);
1008 \endcode
1009
1010 \subsection run_qs_sec Step 4: Run algorithm.
1011 --- starting from the initial iterate @b x, applied to objective function @b obj.
1012
1013 \code
1014 algo.run(x, obj);
1015 \endcode
1016 */
1017
1018 /** @defgroup interface_group User Interface
1019 * \brief ROL's interface, to be implemented by the user.
1020 */
1021
1022 /** @defgroup la_group Linear Algebra Interface
1023 * @ingroup interface_group
1024 * \brief ROL's linear algebra or vector space interface.
1025 */
1026
1027 /** @defgroup func_group Functional Interface
1028 * @ingroup interface_group
1029 * \brief ROL's functional interface.
1030
1031 ROL is used for the numerical solution of smooth optimization problems
1032 \f[
1033 \begin{array}{rl}
1034 \displaystyle \min_{x} & f(x) \\
1035 \mbox{subject to} & c(x) = 0 \,, \\
1036 & a \le x \le b \,,
1037 \end{array}
1038 \f]
1039 where:
1040 \li \f$f : \mathcal{X} \rightarrow \mathbb{R}\f$ is a Fréchet differentiable functional,
1041 \li \f$c : \mathcal{X} \rightarrow \mathcal{C}\f$ is a Fréchet differentiable operator,
1042 \li \f$\mathcal{X}\f$ and \f$\mathcal{C}\f$ are Banach spaces of functions, and
1043 \li \f$a \le x \le b\f$ defines pointwise (componentwise) bounds on \f$x\f$.
1044
1045 This formulation encompasses a variety of useful problem scenarios.
1046
1047 First, the vector spaces \f$\mathcal{X}\f$ and \f$\mathcal{C}\f$, to be defined by the user
1048 through the \ref la_group, can represent real spaces, such as \f$\mathcal{X} = \mathbb{R}^n\f$ and
1049 \f$\mathcal{C} = \mathbb{R}^m\f$, and function spaces, such as Hilbert and Banach function spaces.
1050 ROL's vector space abstractions enable efficent implementations of general optimization algorithms,
1051 spanning traditional nonlinear programming (NLP), optimization with partial differential
1052 equation (PDE) or differential algebraic equation (DAE) constraints, and
1053 stochastic optimization.
1054
1055 Second, ROL's core methods can solve four types of smooth optimization problems, depending on the
1056 structure of the constraints.
1057 \li @b Type-U. No constraints (where \f$c(x) = 0\f$ and \f$a \le x \le b\f$ are absent):
1058 \f[
1059 \begin{array}{rl}
1060 \displaystyle \min_{x} & f(x)
1061 \end{array}
1062 \f]
1063 These problems are known as unconstrained optimization problems.
1064 The user implements the methods of the #ROL::Objective interface.
1065 \li @b Type-B. Bound constraints (where \f$c(x) = 0\f$ is absent):
1066 \f[
1067 \begin{array}{rl}
1068 \displaystyle \min_{x} & f(x) \\
1069 \mbox{subject to} & a \le x \le b \,.
1070 \end{array}
1071 \f]
1072 This problem is typically handled using projections on active sets or primal-dual active-set methods.
1073 ROL provides example implementations of the projections for simple box constraints.
1074 Other projections are defined by the user.
1075 The user implements the methods of the #ROL::BoundConstraint interface.
1076 \li @b Type-E. Equality constraints, generally nonlinear and nonconvex (where \f$a \le x \le b\f$ is absent):
1077 \f[
1078 \begin{array}{rl}
1079 \displaystyle \min_{x} & f(x) \\
1080 \mbox{subject to} & c(x) = 0 \,.
1081 \end{array}
1082 \f]
1083 Equality constraints are handled in ROL using, for example, matrix-free composite-step methods, including sequential quadratic programming (SQP).
1084 The user implements the methods of the #ROL::EqualityConstraint interface.
1085 \li @b Type-EB. Equality and bound constraints:
1086 \f[
1087 \begin{array}{rl}
1088 \displaystyle \min_{x} & f(x) \\
1089 \mbox{subject to} & c(x) = 0 \\
1090 & a \le x \le b \,.
1091 \end{array}
1092 \f]
1093 This formulation includes general inequality constraints.
1094 For example, we can consider the reformulation:
1095 \f[
1096 \begin{array}{rlcccrl}
1097 \displaystyle \min_{x} & f(x) &&&& \displaystyle \min_{x,s} & f(x) \\
1098 \mbox{subject to} & c(x) \le 0 & & \quad \longleftrightarrow \quad & & \mbox{subject to} & c(x) + s = 0 \,, \\
1099 &&&&&& s \ge 0 \,.
1100 \end{array}
1101 \f]
1102 ROL uses a combination of matrix-free composite-step methods, projection methods and primal-dual active set methods to solve these problems.
1103 The user implements the methods of the #ROL::EqualityConstraint and the #ROL::BoundConstraint interfaces.
1104
1105 Third, ROL's design enables streamlined algorithmic extensions, such as the \ref stochastic_group capability.
1106
1107 ---
1108 */
1109
1110 /** @defgroup step_group Optimization Steps
1111 * \brief ROL's optimization steps.
1112 */
1113
1114 /** @defgroup extensions_group Algorithmic Extensions
1115 * \brief ROL's algorithmic extensions.
1116 */
1117
1118 /** @defgroup stochastic_group Stochastic Optimization
1119 * @ingroup extensions_group
1120 * \brief ROL's stochastic optimization capability.
1121 */
1122
1123 /** @defgroup risk_group Risk Measures
1124 * @ingroup stochastic_group
1125 * \brief ROL's risk measure implementations.
1126 */
1127
1128 /** @defgroup dynamic_group Dynamic functions
1129 * @ingroup interface_group
1130 * \brief ROL's interfaces for time-dependent constraints and objectives
1131 */
1132
1133 /** @defgroup examples_group Examples
1134 * \brief ROL's examples
1135 * <ul>
1136 * <li><b>Unconstrained Examples</b>
1137 * <ol>
1138 * <li>\link rol/example/rosenbrock/example_01.cpp Minimizing the Rosenbrock function\endlink</li>
1139 * <li>\link rol/example/zakharov/example_01.cpp Minimizing the Zakharov function\endlink</li>
1140 * <li>\link rol/example/sacado/example_01.hpp Using Sacado with ROL\endlink</li>
1141 * <li>\link rol/example/dual-spaces/rosenbrock-1/example_01.cpp Using Dual Spaces\endlink</li>
1142 * </ol>
1143 * <li><b>Constrained Examples</b></li>
1144 * <ol>
1145 * <li>\link rol/example/dual-spaces/simple-eq-constr-1/example_01.cpp Using Dual Spaces\endlink</li>
1146 * <li>\link rol/example/sacado/example_02.hpp Using Sacado with ROL\endlink</li>
1147 * <li>\link rol/example/poisson-control/example_01.cpp Poisson control\endlink</li>
1148 * <li>\link rol/example/poisson-inversion/example_01.cpp Poisson inversion\endlink</li>
1149 * <li>\link rol/example/burgers-control/example_01.cpp Burgers control\endlink</li>
1150 * <li>\link rol/example/gross-pitaevskii/example_01.hpp Minimizing the Gross-Pitaevskii functional \endlink</li>
1151 * <li>\link rol/example/gross-pitaevskii/example_02.hpp Gross-Pitaevskii functional with \f$H^1\f$ gradient \endlink</li>
1152 * </ol>
1153 * <li><b>Examples using Third Party Libraries</b></li>
1154 * <ol>
1155 * <li>\link rol/example/json/example_01.cpp Using a JSON file to provide ROL with parameters\endlink</li>
1156 * </ol>
1157 * </ul>
1158 */
1159
1160 #endif
1161