13#ifndef OPTIMIST_SOLVER_HH
14#define OPTIMIST_SOLVER_HH
42 template <
typename Real, Integer SolInDim, Integer SolOutDim,
typename DerivedSolver,
bool ForceEigen = false>
47 static_assert(SolInDim >
static_cast<Integer>(0) && SolOutDim >
static_cast<Integer>(0),
48 "Negative-dimensional optimization problem? Are you serious?");
54 using InputType =
typename std::conditional_t<ForceEigen || (SolInDim > 1), Eigen::Vector<Real, SolInDim>, Real>;
55 using OutputType =
typename std::conditional_t<ForceEigen || (SolOutDim > 1), Eigen::Vector<Real, SolOutDim>, Real>;
61 using FirstDerivativeType = std::conditional_t<ForceEigen || (SolInDim > 1) || (SolOutDim > 1), Eigen::Matrix<Real, SolOutDim, SolInDim>, Real>;
63 std::conditional_t<SolInDim == 1 || SolOutDim == 1, Eigen::Matrix<Real, SolInDim, SolInDim>,
64 std::vector<Eigen::Matrix<Real, SolInDim, SolInDim>>>, Real>;
108 if constexpr (ForceEigen || SolInDim > 1) {
109 this->m_lower_bound.setConstant(-INFTY);
110 this->m_upper_bound.setConstant(INFTY);
112 this->m_lower_bound = -INFTY;
113 this->m_upper_bound = INFTY;
115 this->m_trace.reserve(this->m_max_iterations * this->m_max_relaxations);
125 static_cast<const DerivedSolver *
>(
this)->solve_impl(function, x_ini, x_sol);
138 static_cast<const DerivedSolver *
>(
this)->solve_impl(function, first_derivative, x_ini, x_sol);
152 static_cast<const DerivedSolver *
>(
this)->solve_impl(function, first_derivative, second_derivative, x_ini, x_sol);
166 #define CMD "Optimist::Solver::bounds(...): "
168 if constexpr (ForceEigen || SolInDim > 1) {
169 OPTIMIST_ASSERT((this->m_upper_bound - t_lower_bound).minCoeff() <= 0.0,
170 CMD "invalid or degenarate bounds detected.");
173 CMD "invalid or degenarate bounds detected.");
175 this->m_lower_bound = t_lower_bound;
191 #define CMD "Optimist::Solver::bounds(...): "
193 if constexpr (ForceEigen || SolInDim > 1) {
194 OPTIMIST_ASSERT((t_upper_bound - this->m_lower_bound).minCoeff() <= 0.0,
195 CMD "invalid or degenarate bounds detected.");
198 CMD "invalid or degenarate bounds detected.");
200 this->m_upper_bound = t_upper_bound;
212 #define CMD "Optimist::Solver::bounds(...): "
214 if constexpr (ForceEigen || SolInDim > 1) {
216 CMD "invalid or degenarate bounds detected.");
219 CMD "invalid or degenarate bounds detected.");
221 this->m_lower_bound = t_lower_bound;
222 this->m_upper_bound = t_upper_bound;
251 OPTIMIST_ASSERT(!std::isnan(t_max_function_evaluations) && std::isfinite(t_max_function_evaluations),
252 "Optimist::Solver::max_function_evaluations(...): invalid input detected.");
253 this->m_max_function_evaluations = t_max_function_evaluations;
283 "Optimist::Solver::max_first_derivative_evaluations(...): invalid input detected.");
307 "Optimist::Solver::max_second_derivative_evaluations(...): invalid input detected.");
328 OPTIMIST_ASSERT(!std::isnan(t_max_iterations) && std::isfinite(t_max_iterations),
329 "Optimist::Solver::max_iterations(...): invalid input detected.");
330 this->m_max_iterations = t_max_iterations;
337 Real
alpha()
const {
return this->m_alpha;}
346 !std::isnan(t_alpha) && std::isfinite(t_alpha) && 0.0 <= t_alpha && t_alpha <= 1.0,
347 "Optimist::Solver::alpha(...): invalid input detected.");
348 this->m_alpha = t_alpha;
369 OPTIMIST_ASSERT(!std::isnan(t_max_relaxations) && std::isfinite(t_max_relaxations),
370 "Optimist::Solver::max_relaxations(...): invalid input detected.");
371 this->m_max_relaxations = t_max_relaxations;
387 !std::isnan(t_tolerance) && std::isfinite(t_tolerance) && t_tolerance > 0.0,
388 "Optimist::Solver::tolerance(...): invalid input detected.");
389 this->m_tolerance = t_tolerance;
440 std::string
task()
const {
return this->m_task;}
446 void task(std::string t_task) {this->m_task = t_task;}
464 std::ostream &
ostream()
const {
return *this->m_ostream;}
470 void ostream(std::ostream & t_ostream) {this->m_ostream = &t_ostream;}
481 #define CMD "Optimist::Solver::solve(...): "
483 static_assert(DerivedSolver::requires_function,
484 CMD "the solver requires a function.");
485 return static_cast<DerivedSolver *
>(
this)->
solve(function,
nullptr,
nullptr, x_ini, x_sol);
501 #define CMD "Optimist::Solver::solve(...): "
503 static_assert(DerivedSolver::requires_function,
504 CMD "the solver requires a function.");
505 static_assert(DerivedSolver::requires_first_derivative,
506 CMD "the solver requires the first derivative.");
507 return static_cast<DerivedSolver *
>(
this)->
solve(function, first_derivative,
nullptr, x_ini, x_sol);
524 #define CMD "Optimist::Solver::solve(...): "
526 static_assert(DerivedSolver::requires_function,
527 CMD "the solver requires the function.");
528 static_assert(DerivedSolver::requires_first_derivative,
529 CMD "the solver requires the first derivative.");
530 static_assert(DerivedSolver::requires_second_derivative,
531 CMD "the solver requires the second derivative.");
532 return static_cast<DerivedSolver *
>(
this)->
solve(function, first_derivative, second_derivative,
548 template <Integer FunInDim, Integer FunOutDim,
typename DerivedFunction>
552 #define CMD "Optimist::Solver::rootfind(...): "
554 static_assert(SolInDim == FunInDim,
555 CMD "solver input dimension must be equal to the function input dimension.");
556 static_assert(SolOutDim == FunOutDim || SolOutDim == 1,
557 CMD "solver output dimension must be equal to the function output dimension or 1.");
558 static_assert(!(SolInDim == 1 && DerivedSolver::is_optimizer),
559 CMD "one-dimensional optimizers do not support root-finding problems.");
560 return this->
solve(function, x_ini, x_sol, SolOutDim != FunOutDim || (SolInDim > 1 && SolOutDim == 1));
575 template <Integer FunInDim, Integer FunOutDim,
typename DerivedFunction>
579 #define CMD "Optimist::Solver::optimize(...): "
581 static_assert(SolInDim == FunInDim,
582 CMD "solver input dimension must be equal to the function input dimension.");
583 static_assert(SolOutDim == 1,
584 CMD "solver output dimension must be equal to the function output dimension or 1.");
585 static_assert(!(SolInDim == 1 && DerivedSolver::is_rootfinder),
586 CMD "one-dimensional root-finders do not support optimization problems.");
587 return this->
solve(function, x_ini, x_sol,
true);
596 std::string
name()
const {
return static_cast<const DerivedSolver *
>(
this)->name_impl();};
610 template <Integer FunInDim, Integer FunOutDim,
typename DerivedFunction>
614 #define CMD "Optimist::Solver::solve(...): "
618 static_assert(SolInDim == FunInDim,
619 CMD "solver input dimension must be equal to the function input dimension.");
620 static_assert(SolOutDim == FunOutDim || SolOutDim == 1,
621 CMD "solver output dimension must be equal to the function output dimension or 1.");
624 auto function_wrapper = [&function, is_optimization](
const InputType & x,
OutputType & out)
626 typename FunctionType::OutputType f;
628 if (is_optimization) {
629 if constexpr (FunOutDim == 1) {out = 0.5*f*f;}
630 else if constexpr (SolOutDim != FunOutDim) {out = 0.5*f.squaredNorm();}
631 else {
OPTIMIST_ERROR(
CMD "optimization problem with inconsistent output in function.");}
633 if constexpr (SolOutDim == FunOutDim) {
636 else {
OPTIMIST_ERROR(
CMD "root-finding problem with inconsistent output in function.");}
642 typename FunctionType::FirstDerivativeType J; function.
first_derivative(x, J);
644 if (is_optimization) {
645 typename FunctionType::OutputType f; function.
evaluate(x, f);
646 this->m_function_evaluations++;
647 if constexpr (FunInDim == 1 && FunOutDim == 1) {out = J*f;}
648 else if constexpr (SolOutDim != FunOutDim) {out = J.transpose()*f;}
649 else {
OPTIMIST_ERROR(
CMD "optimization problem inconsistent output in first derivative.");}
651 if constexpr (SolOutDim == FunOutDim) {
654 else {
OPTIMIST_ERROR(
CMD "root-finding problem with inconsistent output in first derivative.");}
662 if (is_optimization) {
663 typename FunctionType::OutputType f; function.
evaluate(x, f);
664 typename FunctionType::FirstDerivativeType J; function.
first_derivative(x, J);
665 this->m_function_evaluations++; this->m_first_derivative_evaluations++;
666 if constexpr (FunInDim == 1 && FunOutDim == 1) {out = J*J + f*H;}
667 else if constexpr (SolOutDim != FunOutDim) {
668 out = J.transpose()*J;
669 for (
Integer i = 0; i < static_cast<Integer>(H.size()); ++i) {out += f(i)*H[i];}
671 else {
OPTIMIST_ERROR(
CMD "optimization problem with inconsistent output in second derivative.");}
673 if constexpr (SolOutDim == FunOutDim) {out = H;}
674 else {
OPTIMIST_ERROR(
CMD "root-finding problem with inconsistent output in second derivative.");}
679 if constexpr (DerivedSolver::requires_function && !DerivedSolver::requires_first_derivative &&
680 !DerivedSolver::requires_second_derivative) {
681 return static_cast<DerivedSolver *
>(
this)->
solve(function_wrapper, x_ini, x_sol);
682 }
else if constexpr (DerivedSolver::requires_function && DerivedSolver::requires_first_derivative &&
683 !DerivedSolver::requires_second_derivative) {
684 return static_cast<DerivedSolver *
>(
this)->
solve(function_wrapper, first_derivative_wrapper,
686 }
else if constexpr (DerivedSolver::requires_function && DerivedSolver::requires_first_derivative &&
687 DerivedSolver::requires_second_derivative) {
688 return static_cast<DerivedSolver *
>(
this)->
solve(function_wrapper, first_derivative_wrapper,
689 second_derivative_wrapper, x_ini, x_sol);
702 this->m_function_evaluations =
static_cast<Integer>(0);
703 this->m_first_derivative_evaluations =
static_cast<Integer>(0);
704 this->m_second_derivative_evaluations =
static_cast<Integer>(0);
705 this->m_iterations =
static_cast<Integer>(0);
706 this->m_relaxations =
static_cast<Integer>(0);
707 this->m_converged =
false;
708 this->m_trace.clear();
720 "Optimist::" << this->
name() <<
"::evaluate_function(...): maximum allowed function evaluations reached.");
721 ++this->m_function_evaluations;
734 "Optimist::" << this->
name() <<
"::evaluate_first_derivative(...): maximum allowed first derivative evaluations reached.");
735 ++this->m_first_derivative_evaluations;
748 "Optimist::" << this->
name() <<
"::evaluate_second_derivative(...): maximum allowed second derivative evaluations reached.");
749 ++this->m_second_derivative_evaluations;
773 Real step_norm_old, step_norm_new, residuals_old, residuals_new, tau{1.0};
774 for (this->m_relaxations =
static_cast<Integer>(0); this->m_relaxations < this->m_max_relaxations; ++this->m_relaxations)
777 step_new = tau * step_old;
778 x_new = x_old + step_new;
782 if constexpr (ForceEigen || (SolInDim > 1 && SolOutDim > 1)) {
783 residuals_old = function_old.norm();
784 residuals_new = function_new.norm();
785 step_norm_old = step_old.norm();
786 step_norm_new = step_new.norm();
788 residuals_old = std::abs(function_old);
789 residuals_new = std::abs(function_new);
790 step_norm_old = std::abs(step_old);
791 step_norm_new = std::abs(step_new);
793 if (residuals_new < residuals_old || step_norm_new < (Real(1.0)-tau/Real(2.0))*step_norm_old) {
796 tau *= this->m_alpha;
823 << c_tl << h_78 << c_tr << std::endl
824 << v_ll <<
"Solver:" << std::setw(69) << this->
name() << v_rr << std::endl
825 << v_ll <<
"Task: " << std::setw(69) << this->
task() << v_rr << std::endl
826 << j_ll << h_7 << j_tt << h_7 << j_tt << h_7 << j_tt << h_7 << j_tt << h_7 << j_tt << h_14 << j_tt << h_23 << j_rr << std::endl
827 << v_ll <<
"#Iter" << v_lc <<
" #f" << v_lc <<
" #Df" << v_lc <<
" #DDf" << v_lc <<
" #Rlx" << v_lc
828 <<
" ║f(x)║₂ " << v_lc <<
"Additional notes" << std::setw(9) << v_rr << std::endl
829 << j_ll << h_7 << j_cc << h_7 << j_cc << h_7 << j_cc << h_7 << j_cc << h_7 << j_cc << h_14 << j_cc << h_23 << j_rr << std::endl;
851 << j_ll << h_7 << j_bb << h_7 << j_bb << h_7 << j_bb << h_7 << j_bb << h_7 << j_bb << h_14 << j_bb << h_23 << j_rr << std::endl
852 << v_ll << std::setw(40) << (this->m_converged ?
"CONVERGED" :
"NOT CONVERGED") << std::setw(40) << v_rr << std::endl
853 << c_bl << h_78 << c_br << std::endl;
860 void info(Real residuals, std::string
const & notes =
"-")
866 *this->m_ostream << v_ll
867 << std::setw(5) << this->m_iterations << v_lc
868 << std::setw(5) << this->m_function_evaluations << v_lc
869 << std::setw(5) << this->m_first_derivative_evaluations << v_lc
870 << std::setw(5) << this->m_second_derivative_evaluations << v_lc
871 << std::setw(5) << this->m_relaxations << v_lc
872 << std::setw(12) << std::scientific << std::setprecision(6) << residuals << v_lc
873 << notes << std::setw(25-notes.length()) << v_rr << std::endl;
#define OPTIMIST_ERROR(MSG)
Definition Optimist.hh:33
#define OPTIMIST_BASIC_CONSTANTS(Real)
Definition Optimist.hh:70
#define OPTIMIST_ASSERT(COND, MSG)
Definition Optimist.hh:43
Class container for the generic function.
Definition Function.hh:43
void second_derivative(const InputType &x, SecondDerivativeType &out) const
Definition Function.hh:102
void first_derivative(const InputType &x, FirstDerivativeType &out) const
Definition Function.hh:92
void evaluate(const InputType &x, OutputType &out) const
Definition Function.hh:82
std::conditional_t< ForceEigen||(SolInDim > 1)||(SolOutDim > 1), Eigen::Matrix< Real, SolOutDim, SolInDim >, Real > FirstDerivativeType
Definition SolverBase.hh:61
Integer max_iterations() const
Definition SolverBase.hh:321
typename std::vector< InputType > TraceType
Definition SolverBase.hh:58
bool m_damped
Definition SolverBase.hh:95
void max_second_derivative_evaluations(Integer second_derivative_evaluations)
Definition SolverBase.hh:303
void ostream(std::ostream &t_ostream)
Definition SolverBase.hh:470
bool converged() const
Definition SolverBase.hh:452
bool rootfind(FunctionBase< Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen &&FunOutDim==1 > const &function, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:549
const InputType & upper_bound() const
Definition SolverBase.hh:184
void alpha(Real t_alpha)
Definition SolverBase.hh:343
Integer max_first_derivative_evaluations() const
Definition SolverBase.hh:273
void max_iterations(Integer t_max_iterations)
Definition SolverBase.hh:327
void info(Real residuals, std::string const ¬es="-")
Definition SolverBase.hh:860
void lower_bound(const InputType &t_lower_bound)
Definition SolverBase.hh:165
constexpr Integer output_dimension() const
Definition SolverBase.hh:237
const TraceType & trace() const
Definition SolverBase.hh:458
Integer m_iterations
Definition SolverBase.hh:86
void bounds(const InputType &t_lower_bound, const InputType &t_upper_bound)
Definition SolverBase.hh:210
Integer max_relaxations() const
Definition SolverBase.hh:361
Real alpha() const
Definition SolverBase.hh:337
bool solve(FunctionWrapper function, FirstDerivativeWrapper first_derivative, SecondDerivativeWrapper second_derivative, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:521
SolverBase(FunctionWrapper function, FirstDerivativeWrapper first_derivative, SecondDerivativeWrapper second_derivative, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:149
Integer m_max_iterations
Definition SolverBase.hh:87
std::string task() const
Definition SolverBase.hh:440
SolverBase(FunctionWrapper function, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:124
void enable_damped_mode()
Definition SolverBase.hh:429
void upper_bound(const InputType &t_upper_bound)
Definition SolverBase.hh:190
Integer relaxations() const
Definition SolverBase.hh:355
typename std::function< void(const InputType &, SecondDerivativeType &)> SecondDerivativeWrapper
Definition SolverBase.hh:69
void verbose_mode(bool t_verbose)
Definition SolverBase.hh:396
SolverBase(FunctionWrapper function, FirstDerivativeWrapper first_derivative, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:135
std::string name() const
Definition SolverBase.hh:596
void header()
Definition SolverBase.hh:806
Integer max_second_derivative_evaluations() const
Definition SolverBase.hh:297
bool solve(FunctionWrapper function, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:479
bool optimize(FunctionBase< Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen &&FunOutDim==1 > const &function, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:576
Integer function_evaluations() const
Definition SolverBase.hh:243
Real tolerance() const
Definition SolverBase.hh:378
Integer iterations() const
Definition SolverBase.hh:315
typename std::conditional_t< ForceEigen||(SolInDim > 1), Eigen::Vector< Real, SolInDim >, Real > InputType
Definition SolverBase.hh:54
void reset()
Definition SolverBase.hh:700
void max_relaxations(Integer t_max_relaxations)
Definition SolverBase.hh:367
std::ostream & ostream() const
Definition SolverBase.hh:464
void damped_mode(bool t_damped)
Definition SolverBase.hh:418
Integer m_max_first_derivative_evaluations
Definition SolverBase.hh:82
void evaluate_second_derivative(SecondDerivativeWrapper function, const InputType &x, SecondDerivativeType &out)
Definition SolverBase.hh:745
bool m_verbose
Definition SolverBase.hh:94
typename std::function< void(const InputType &, OutputType &)> FunctionWrapper
Definition SolverBase.hh:67
Integer m_max_second_derivative_evaluations
Definition SolverBase.hh:83
void evaluate_first_derivative(FirstDerivativeWrapper function, const InputType &x, FirstDerivativeType &out)
Definition SolverBase.hh:731
bool m_converged
Definition SolverBase.hh:100
Real m_alpha
Definition SolverBase.hh:88
InputType m_upper_bound
Definition SolverBase.hh:73
void enable_verbose_mode()
Definition SolverBase.hh:407
void disable_damped_mode()
Definition SolverBase.hh:434
constexpr Integer input_dimension() const
Definition SolverBase.hh:231
Integer max_function_evaluations() const
Definition SolverBase.hh:260
typename std::conditional_t< ForceEigen||(SolOutDim > 1), Eigen::Vector< Real, SolOutDim >, Real > OutputType
Definition SolverBase.hh:55
TraceType m_trace
Definition SolverBase.hh:101
bool verbose_mode() const
Definition SolverBase.hh:402
const InputType & lower_bound() const
Definition SolverBase.hh:159
bool damp(FunctionWrapper function, InputType const &x_old, InputType const &function_old, InputType const &step_old, InputType &x_new, InputType &function_new, InputType &step_new)
Definition SolverBase.hh:770
Integer second_derivative_evaluations() const
Definition SolverBase.hh:291
void store_trace(const InputType &x)
Definition SolverBase.hh:757
void bottom()
Definition SolverBase.hh:836
SolverBase()
Definition SolverBase.hh:107
std::ostream * m_ostream
Definition SolverBase.hh:96
Integer m_function_evaluations
Definition SolverBase.hh:76
bool solve(FunctionBase< Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen &&FunOutDim==1 > const &function, const InputType &x_ini, InputType &x_sol, bool is_optimization)
Definition SolverBase.hh:611
std::conditional_t< ForceEigen||(SolInDim > 1)||(SolOutDim > 1), std::conditional_t< SolInDim==1||SolOutDim==1, Eigen::Matrix< Real, SolInDim, SolInDim >, std::vector< Eigen::Matrix< Real, SolInDim, SolInDim > > >, Real > SecondDerivativeType
Definition SolverBase.hh:62
Integer m_first_derivative_evaluations
Definition SolverBase.hh:77
Integer m_relaxations
Definition SolverBase.hh:89
void max_first_derivative_evaluations(Integer first_derivative_evaluations)
Definition SolverBase.hh:279
Real m_tolerance
Definition SolverBase.hh:93
bool solve(FunctionWrapper function, FirstDerivativeWrapper first_derivative, const InputType &x_ini, InputType &x_sol)
Definition SolverBase.hh:498
std::string m_task
Definition SolverBase.hh:99
bool damped_mode() const
Definition SolverBase.hh:424
Integer m_second_derivative_evaluations
Definition SolverBase.hh:78
Integer m_max_relaxations
Definition SolverBase.hh:90
Integer first_derivative_evaluations() const
Definition SolverBase.hh:267
void disable_verbose_mode()
Definition SolverBase.hh:412
void evaluate_function(FunctionWrapper function, const InputType &x, OutputType &out)
Definition SolverBase.hh:717
void tolerance(Real t_tolerance)
Definition SolverBase.hh:385
Integer m_max_function_evaluations
Definition SolverBase.hh:81
void max_function_evaluations(Integer t_max_function_evaluations)
Definition SolverBase.hh:249
void task(std::string t_task)
Definition SolverBase.hh:446
typename std::function< void(const InputType &, FirstDerivativeType &)> FirstDerivativeWrapper
Definition SolverBase.hh:68
InputType m_lower_bound
Definition SolverBase.hh:72
Namespace for the Optimist library.
Definition Optimist.hh:87
static std::string table_vertical_line()
Retrieve the Unicode character for the vertical line of a table.
Definition Optimist.hh:173
static std::string table_bottom_right_corner()
Retrieve the Unicode character for the bottom-right corner of a table.
Definition Optimist.hh:119
static std::string table_top_left_corner()
Retrieve the Unicode character for the top-left corner of a table.
Definition Optimist.hh:101
static std::string table_bottom_left_corner()
Retrieve the Unicode character for the bottom-left corner of a table.
Definition Optimist.hh:113
OPTIMIST_DEFAULT_INTEGER_TYPE Integer
The Integer type as used for the API.
Definition Optimist.hh:95
static std::string table_center_cross()
Retrieve the Unicode character for the center cross of a table.
Definition Optimist.hh:149
static std::string table_bottom_junction()
Retrieve the Unicode character for the bottom junction of a table.
Definition Optimist.hh:143
static std::string table_top_junction()
Retrieve the Unicode character for the top junction of a table.
Definition Optimist.hh:137
static std::string table_top_right_corner()
Retrieve the Unicode character for the top-right corner of a table.
Definition Optimist.hh:107
static std::string table_horizontal_line()
Retrieve the Unicode character for the horizontal line of a table.
Definition Optimist.hh:155
static std::string table_right_junction()
Retrieve the Unicode character for the right junction of a table.
Definition Optimist.hh:131
static std::string table_left_junction()
Retrieve the Unicode character for the left junction of a table.
Definition Optimist.hh:125