SolverBase.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco, Mattia Piazza and Enrico Bertolazzi.                       *
 *                                                                                               *
 * The Optimist project is distributed under the BSD 2-Clause License.                           *
 *                                                                                               *
 * Davide Stocco                          Mattia Piazza                        Enrico Bertolazzi *
 * University of Trento               University of Trento                  University of Trento *
 * davide.stocco@unitn.it            mattia.piazza@unitn.it           enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef OPTIMIST_SOLVER_HH
#define OPTIMIST_SOLVER_HH
 
#include "Optimist.hh"
#include "Optimist/Function.hh"
 
namespace Optimist
{
 
  /*\
   |   ____        _                ____
   |  / ___|  ___ | |_   _____ _ __| __ )  __ _ ___  ___
   |  \___ \ / _ \| \ \ / / _ \ '__|  _ \ / _` / __|/ _ \
   |   ___) | (_) | |\ V /  __/ |  | |_) | (_| \__ \  __/
   |  |____/ \___/|_| \_/ \___|_|  |____/ \__,_|___/\___|
   |
  \*/

  template <typename Real, Integer SolInDim, Integer SolOutDim, typename DerivedSolver, bool ForceEigen = false>
  class SolverBase
  {
  public:
    // Fancy static assertions (just for fun, don't take it too seriously)
    static_assert(SolInDim > 0 && SolOutDim > 0,
      "Negative-dimensional optimization problem? Are you serious?");
 
    OPTIMIST_BASIC_CONSTANTS(Real)
 
    // I/O types
    using InputType  = typename std::conditional_t<ForceEigen || (SolInDim > 1),
      Eigen::Vector<Real, SolInDim>, Real>;
    using OutputType = typename std::conditional_t<ForceEigen || (SolOutDim > 1),
      Eigen::Vector<Real, SolOutDim>, Real>;
 
    // Trace types
    using TraceType = typename std::vector<InputType>;
 
    // Derivative types
    using FirstDerivativeType = std::conditional_t<ForceEigen || (SolInDim > 1) || (SolOutDim > 1),
      Eigen::Matrix<Real, SolOutDim, SolInDim>, Real>;
    using SecondDerivativeType = 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>;
 
  protected:
    // Bounds (may not be used)
    InputType m_lower_bound; 
    InputType m_upper_bound; 
 
    // Evaluations
    Integer m_function_evaluations{0};          
    Integer m_first_derivative_evaluations{0};  
    Integer m_second_derivative_evaluations{0}; 
 
    // Maximum allowed evaluations
    Integer m_max_function_evaluations{1000};          
    Integer m_max_first_derivative_evaluations{1000};  
    Integer m_max_second_derivative_evaluations{1000}; 
 
    // Iterations and relaxations
    Integer m_iterations{0};       
    Integer m_max_iterations{100}; 
    Real    m_alpha{0.8};          
    Integer m_relaxations{0};      
    Integer m_max_relaxations{10}; 
 
    // Settings
    Real m_tolerance{EPSILON_LOW};        
    bool m_verbose{false};                
    bool m_damped{true};                  
    std::ostream * m_ostream{&std::cout}; 
 
    // Convergence output flag and trace
    std::string m_task{"Undefined"}; 
    bool        m_converged{false};  
    TraceType   m_trace;             
 
  public:
    SolverBase() {
      if constexpr (ForceEigen || SolInDim > 1) {
        this->m_lower_bound.setConstant(-INFTY);
        this->m_upper_bound.setConstant(INFTY);
      } else {
        this->m_lower_bound = -INFTY;
        this->m_upper_bound = INFTY;
      }
      this->m_trace.reserve(this->m_max_iterations * this->m_max_relaxations);
    }

    template <typename FunctionLambda>
    SolverBase(FunctionLambda && function, InputType const & x_ini, InputType & x_sol) : SolverBase() {
      static_cast<const DerivedSolver *>(this)->solve_impl(
        std::forward<FunctionLambda>(function),
        x_ini, x_sol);
    }

    template <typename FunctionLambda, typename FirstDerivativeLambda>
    SolverBase(FunctionLambda && function, FirstDerivativeLambda && first_derivative, InputType const & x_ini,
      InputType & x_sol) : SolverBase()
    {
      static_cast<const DerivedSolver *>(this)->solve_impl(
        std::forward<FunctionLambda>(function),
        std::forward<FirstDerivativeLambda>(first_derivative),
        x_ini, x_sol);
    }

    template <typename FunctionLambda, typename FirstDerivativeLambda, typename SecondDerivativeLambda>
    SolverBase(FunctionLambda && function, FirstDerivativeLambda && first_derivative, SecondDerivativeLambda
      && second_derivative, InputType const & x_ini, InputType & x_sol) : SolverBase()
    {
      static_cast<const DerivedSolver *>(this)->solve_impl(
        std::forward<FunctionLambda>(function),
        std::forward<FirstDerivativeLambda>(first_derivative),
        std::forward<SecondDerivativeLambda>(second_derivative),
        x_ini, x_sol);
    }

    InputType const & lower_bound() const {return this->m_lower_bound;}

    void lower_bound(InputType const & t_lower_bound) {
      #define CMD "Optimist::Solver::bounds(...): "
 
      if constexpr (ForceEigen || SolInDim > 1) {
        OPTIMIST_ASSERT((this->m_upper_bound - t_lower_bound).minCoeff() <= 0.0,
          CMD "invalid or degenarate bounds detected.");
      } else {
        OPTIMIST_ASSERT(this->m_upper_bound > t_lower_bound,
          CMD "invalid or degenarate bounds detected.");
      }
      this->m_lower_bound = t_lower_bound;
 
      #undef CMD
    }

    InputType const & upper_bound() const {return this->m_upper_bound;}

    void upper_bound(InputType const & t_upper_bound) {
      #define CMD "Optimist::Solver::bounds(...): "
 
      if constexpr (ForceEigen || SolInDim > 1) {
        OPTIMIST_ASSERT((t_upper_bound - this->m_lower_bound).minCoeff() <= 0.0,
          CMD "invalid or degenarate bounds detected.");
      } else {
        OPTIMIST_ASSERT(t_upper_bound > this->m_lower_bound,
          CMD "invalid or degenarate bounds detected.");
      }
      this->m_upper_bound = t_upper_bound;
 
      #undef CMD
    }

    void bounds(InputType const & t_lower_bound, InputType const & t_upper_bound)
    {
      #define CMD "Optimist::Solver::bounds(...): "
 
      if constexpr (ForceEigen || SolInDim > 1) {
        OPTIMIST_ASSERT((t_upper_bound - t_lower_bound).minCoeff() <= 0.0,
          CMD "invalid or degenarate bounds detected.");
      } else {
        OPTIMIST_ASSERT(t_upper_bound > t_lower_bound,
          CMD "invalid or degenarate bounds detected.");
      }
      this->m_lower_bound = t_lower_bound;
      this->m_upper_bound = t_upper_bound;
 
      #undef CMD
    }

    constexpr Integer input_dimension() const {return SolInDim;}

    constexpr Integer output_dimension() const {return SolOutDim;}

    Integer function_evaluations() const {return this->m_function_evaluations;}

    void max_function_evaluations(Integer t_max_function_evaluations)
    {
      OPTIMIST_ASSERT(t_max_function_evaluations > 0,
        "Optimist::Solver::max_function_evaluations(...): invalid input detected.");
      this->m_max_function_evaluations = t_max_function_evaluations;
    }

    Integer max_function_evaluations() const {return this->m_max_function_evaluations;}
 
  protected:
    Integer first_derivative_evaluations() const {return this->m_first_derivative_evaluations;}

    Integer max_first_derivative_evaluations() const {return this->m_max_first_derivative_evaluations;}

    void max_first_derivative_evaluations(Integer first_derivative_evaluations)
    {
      OPTIMIST_ASSERT(first_derivative_evaluations > 0,
        "Optimist::Solver::max_first_derivative_evaluations(...): invalid input detected.");
      this->m_max_first_derivative_evaluations = first_derivative_evaluations;
    }

    Integer second_derivative_evaluations() const {return this->m_second_derivative_evaluations;}

    Integer max_second_derivative_evaluations() const {return this->m_max_second_derivative_evaluations;}

    void max_second_derivative_evaluations(Integer second_derivative_evaluations)
    {
      OPTIMIST_ASSERT(second_derivative_evaluations > 0,
        "Optimist::Solver::max_second_derivative_evaluations(...): invalid input detected.");
      this->m_max_second_derivative_evaluations = second_derivative_evaluations;
    }
 
  public:
    Integer iterations() const {return this->m_iterations;}

    Integer max_iterations() const {return this->max_iterations;}

    void max_iterations(Integer t_max_iterations) {
      OPTIMIST_ASSERT(t_max_iterations > 0,
        "Optimist::Solver::max_iterations(...): invalid input detected.");
      this->m_max_iterations = t_max_iterations;
    }

    Real alpha() const {return this->m_alpha;}

    void alpha(Real t_alpha)
    {
      OPTIMIST_ASSERT(
        !std::isnan(t_alpha) && std::isfinite(t_alpha) && 0.0 <= t_alpha && t_alpha <= 1.0,
        "Optimist::Solver::alpha(...): invalid input detected.");
      this->m_alpha = t_alpha;
    }

    Integer relaxations() const {return this->m_relaxations;}

    Integer max_relaxations() const {return this->max_relaxations;}

    void max_relaxations(Integer t_max_relaxations)
    {
      OPTIMIST_ASSERT(t_max_relaxations > 0,
        "Optimist::Solver::max_relaxations(...): invalid input detected.");
      this->m_max_relaxations = t_max_relaxations;
    }

    Real tolerance() const {return this->m_tolerance;}

    void tolerance(Real t_tolerance) {
      OPTIMIST_ASSERT(!std::isnan(t_tolerance) && std::isfinite(t_tolerance) && t_tolerance > 0.0,
        "Optimist::Solver::tolerance(...): invalid input detected.");
      this->m_tolerance = t_tolerance;
    }

    void verbose_mode(bool t_verbose) {this->m_verbose = t_verbose;}

    bool verbose_mode() const {return this->m_verbose;}

    void enable_verbose_mode() {this->m_verbose = true;}

    void disable_verbose_mode() {this->m_verbose = false;}

    void damped_mode(bool t_damped) {this->m_damped = t_damped;}

    bool damped_mode() const {return this->m_damped;}

    void enable_damped_mode() {this->m_damped = true;}

    void disable_damped_mode() {this->m_damped = false;}

    std::string task() const {return this->m_task;}

    void task(std::string t_task) {this->m_task = t_task;}

    bool converged() const {return this->m_converged;}

    const TraceType & trace() const {return this->m_trace;}

    std::ostream & ostream() const {return *this->m_ostream;}

    void ostream(std::ostream & t_ostream) {this->m_ostream = &t_ostream;}

    template <typename FunctionLambda>
    bool solve(FunctionLambda && function, InputType const & x_ini, InputType & x_sol)
    {
      #define CMD "Optimist::Solver::solve(...): "
 
      static_assert(DerivedSolver::requires_function,
        CMD "the solver requires a function.");
      return static_cast<DerivedSolver *>(this)->solve(
        std::forward<FunctionLambda>(function),
        nullptr, nullptr, x_ini, x_sol);
 
      #undef CMD
    }

    template <typename FunctionLambda, typename FirstDerivativeLambda>
    bool solve(FunctionLambda && function, FirstDerivativeLambda && first_derivative, InputType const & x_ini,
      InputType & x_sol)
    {
      #define CMD "Optimist::Solver::solve(...): "
 
      static_assert(DerivedSolver::requires_function,
        CMD "the solver requires a function.");
      static_assert(DerivedSolver::requires_first_derivative,
        CMD "the solver requires the first derivative.");
      return static_cast<DerivedSolver *>(this)->solve(
        std::forward<FunctionLambda>(function),
        std::forward<FirstDerivativeLambda>(first_derivative),
        nullptr, x_ini, x_sol);
 
      #undef CMD
    }

    template <typename FunctionLambda, typename FirstDerivativeLambda, typename SecondDerivativeLambda>
    bool solve(FunctionLambda && function, FirstDerivativeLambda && first_derivative, SecondDerivativeLambda
      && second_derivative, InputType const & x_ini, InputType & x_sol)
      {
      #define CMD "Optimist::Solver::solve(...): "
 
      static_assert(DerivedSolver::requires_function,
        CMD "the solver requires the function.");
      static_assert(DerivedSolver::requires_first_derivative,
        CMD "the solver requires the first derivative.");
      static_assert(DerivedSolver::requires_second_derivative,
        CMD "the solver requires the second derivative.");
      return static_cast<DerivedSolver *>(this)->solve(
        std::forward<FunctionLambda>(function),
        std::forward<FirstDerivativeLambda>(first_derivative),
        std::forward<SecondDerivativeLambda>(second_derivative),
        x_ini, x_sol);
 
      #undef CMD
    }

    template <Integer FunInDim, Integer FunOutDim, typename DerivedFunction>
    bool rootfind(FunctionBase<Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen && FunOutDim == 1> const & function,
      InputType const & x_ini, InputType & x_sol)
    {
      #define CMD "Optimist::Solver::rootfind(...): "
 
      static_assert(SolInDim == FunInDim,
        CMD "solver input dimension must be equal to the function input dimension.");
      static_assert(SolOutDim == FunOutDim || SolOutDim == 1,
        CMD "solver output dimension must be equal to the function output dimension or 1.");
      static_assert(!(SolInDim == 1 && DerivedSolver::is_optimizer),
        CMD "one-dimensional optimizers do not support root-finding problems.");
      return this->solve(function, x_ini, x_sol, SolOutDim != FunOutDim || (SolInDim > 1 && SolOutDim == 1));
 
      #undef CMD
    }

    template <Integer FunInDim, Integer FunOutDim, typename DerivedFunction>
    bool optimize(FunctionBase<Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen && FunOutDim == 1> const & function,
      InputType const & x_ini, InputType & x_sol)
    {
      #define CMD "Optimist::Solver::optimize(...): "
 
      static_assert(SolInDim == FunInDim,
        CMD "solver input dimension must be equal to the function input dimension.");
      static_assert(SolOutDim == 1,
        CMD "solver output dimension must be equal to the function output dimension or 1.");
      static_assert(!(SolInDim == 1 && DerivedSolver::is_rootfinder),
        CMD "one-dimensional root-finders do not support optimization problems.");
      return this->solve(function, x_ini, x_sol, true);
 
      #undef CMD
    }

    std::string name() const {return static_cast<const DerivedSolver *>(this)->name_impl();};
 
  protected:
    template <Integer FunInDim, Integer FunOutDim, typename DerivedFunction>
    bool solve(FunctionBase<Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen && FunOutDim == 1>
      const & function, InputType const & x_ini, InputType & x_sol, bool is_optimization)
    {
      #define CMD "Optimist::Solver::solve(...): "
 
      using FunctionType = FunctionBase<Real, FunInDim, FunOutDim, DerivedFunction, ForceEigen && FunOutDim == 1>;
 
      static_assert(SolInDim == FunInDim,
        CMD "solver input dimension must be equal to the function input dimension.");
      static_assert(SolOutDim == FunOutDim || SolOutDim == 1,
        CMD "solver output dimension must be equal to the function output dimension or 1.");
 
      // Lambda generators for function and derivatives
      auto function_lambda = [&function, is_optimization] (InputType const & x, OutputType & out) -> bool
      {
        bool success{false};
        typename FunctionType::OutputType f; success = function.evaluate(x, f);
        OPTIMIST_ASSERT(success,
          CMD "function evaluation failed during function computation.");
 
        if (is_optimization) {
          if constexpr (FunOutDim == 1) {out = 0.5*f*f;}
          else if constexpr (SolOutDim != FunOutDim) {out = 0.5*f.squaredNorm();}
          else {OPTIMIST_ERROR(CMD "optimization problem with inconsistent output in function.");}
        } else {
          if constexpr (SolOutDim == FunOutDim) {out = f;}
          else {OPTIMIST_ERROR(CMD "root-finding problem with inconsistent output in function.");}
        }
        return success;
      };
 
      auto first_derivative_lambda = [&function, is_optimization, this] (InputType const & x,
        FirstDerivativeType & out) -> bool
      {
        bool success{false};
        typename FunctionType::FirstDerivativeType J; success = function.first_derivative(x, J);
        OPTIMIST_ASSERT(success,
          CMD "first derivative evaluation failed during first derivative computation.");
 
        if (is_optimization) {
          bool success{false};
          typename FunctionType::OutputType f; success = function.evaluate(x, f);
          OPTIMIST_ASSERT(success,
            CMD "function evaluation failed during first derivative computation.");
          this->m_function_evaluations++;
          if constexpr (FunInDim == 1 && FunOutDim == 1) {out = J*f;}
          else if constexpr (SolOutDim != FunOutDim) {out = J.transpose()*f;}
          else {OPTIMIST_ERROR(CMD "optimization problem inconsistent output in first derivative.");}
        } else {
          if constexpr (SolOutDim == FunOutDim) {out = J;}
          else {OPTIMIST_ERROR(CMD "root-finding problem with inconsistent output in first derivative.");}
        }
        return success;
      };
 
      auto second_derivative_lambda = [&function, is_optimization, this] (InputType const & x,
        SecondDerivativeType & out) -> bool
      {
        bool success{false};
        typename FunctionType::SecondDerivativeType H; success = function.second_derivative(x, H);
        OPTIMIST_ASSERT(success,
          CMD "function evaluation failed during second derivative computation.");
 
        if (is_optimization) {
          typename FunctionType::OutputType f; success = function.evaluate(x, f);
          this->m_function_evaluations++;
          OPTIMIST_ASSERT(success,
            CMD "function evaluation failed during second derivative computation.");
          typename FunctionType::FirstDerivativeType J; success = function.first_derivative(x, J);
          this->m_first_derivative_evaluations++;
          OPTIMIST_ASSERT(success,
            CMD "first derivative evaluation failed during second derivative computation.");
          if constexpr (FunInDim == 1 && FunOutDim == 1) {out = J*J + f*H;}
          else if constexpr (SolOutDim != FunOutDim) {
            out = J.transpose()*J;
            for (Integer i{0}; i < static_cast<Integer>(H.size()); ++i) {out += f(i)*H[i];}
          }
          else {OPTIMIST_ERROR(CMD "optimization problem with inconsistent output in second derivative.");}
        } else {
          if constexpr (SolOutDim == FunOutDim) {out = H;}
          else {OPTIMIST_ERROR(CMD "root-finding problem with inconsistent output in second derivative.");}
        }
        return success;
      };
 
      // Select solver method based on derivative requirements
      if constexpr (DerivedSolver::requires_function && !DerivedSolver::requires_first_derivative &&
        !DerivedSolver::requires_second_derivative) {
        return static_cast<DerivedSolver *>(this)->solve(function_lambda, x_ini, x_sol);
      } else if constexpr (DerivedSolver::requires_function && DerivedSolver::requires_first_derivative &&
        !DerivedSolver::requires_second_derivative) {
        return static_cast<DerivedSolver *>(this)->solve(function_lambda, first_derivative_lambda,
          x_ini, x_sol);
      } else if constexpr (DerivedSolver::requires_function && DerivedSolver::requires_first_derivative &&
        DerivedSolver::requires_second_derivative) {
        return static_cast<DerivedSolver *>(this)->solve(function_lambda, first_derivative_lambda,
          second_derivative_lambda, x_ini, x_sol);
      } else {
        OPTIMIST_ERROR(CMD "no matching function signature found for the solver.");
      }
 
      #undef CMD
    }

    void reset()
    {
      this->m_function_evaluations          = 0;
      this->m_first_derivative_evaluations  = 0;
      this->m_second_derivative_evaluations = 0;
      this->m_iterations                    = 0;
      this->m_relaxations                   = 0;
      this->m_converged                     = false;
      this->m_trace.clear();
    }

    template <typename FunctionLambda>
    bool evaluate_function(FunctionLambda && function, InputType const & x, OutputType & out)
    {
      OPTIMIST_ASSERT(this->m_function_evaluations < this->m_max_function_evaluations,
        "Optimist::" << this-> name() << "::evaluate_function(...): maximum allowed function evaluations reached.");
      ++this->m_function_evaluations;
      return function(x, out);
    }

    template <typename FirstDerivativeLambda>
    bool evaluate_first_derivative(FirstDerivativeLambda && function, InputType const & x, FirstDerivativeType & out)
    {
      OPTIMIST_ASSERT(this->m_first_derivative_evaluations < this->m_max_first_derivative_evaluations,
        "Optimist::" << this-> name() << "::evaluate_first_derivative(...): maximum allowed first derivative evaluations reached.");
      ++this->m_first_derivative_evaluations;
      return function(x, out);
    }

    template <typename SecondDerivativeLambda>
    bool evaluate_second_derivative(SecondDerivativeLambda && function, InputType const & x, SecondDerivativeType & out)
    {
      OPTIMIST_ASSERT(this->m_second_derivative_evaluations < this->m_max_second_derivative_evaluations,
        "Optimist::" << this-> name() << "::evaluate_second_derivative(...): maximum allowed second derivative evaluations reached.");
      ++this->m_second_derivative_evaluations;
      return function(x, out);
    }

    void store_trace(InputType const & x) {this->m_trace.push_back(x);}

    template <typename FunctionLambda>
    bool damp(FunctionLambda && function, InputType const & x_old, InputType const & function_old,
      InputType const & step_old, InputType & x_new, InputType & function_new, InputType & step_new)
    {
      #define CMD "Optimist::Solver::damp(...): "
 
      Real step_norm_old, step_norm_new, residuals_old, residuals_new, tau{1.0};
      for (this->m_relaxations = 0; this->m_relaxations < this->m_max_relaxations; ++this->m_relaxations)
      {
        // Update point
        step_new = tau * step_old;
        x_new = x_old + step_new;
        bool success{this->evaluate_function(std::forward<FunctionLambda>(function), x_new, function_new)};
        OPTIMIST_ASSERT(success,
          CMD "function evaluation failed during damping.");
 
        // Check relaxation
        if constexpr (ForceEigen || (SolInDim > 1 && SolOutDim > 1)) {
          residuals_old = function_old.norm();
          residuals_new = function_new.norm();
          step_norm_old = step_old.norm();
          step_norm_new = step_new.norm();
        } else {
          residuals_old = std::abs(function_old);
          residuals_new = std::abs(function_new);
          step_norm_old = std::abs(step_old);
          step_norm_new = std::abs(step_new);
        }
        if (residuals_new < residuals_old || step_norm_new < (Real(1.0)-tau/Real(2.0))*step_norm_old) {
          return true;
        } else {
          tau *= this->m_alpha;
        }
      }
      return false;
 
      #undef CMD
    }

    void header()
    {
      std::string c_tl{table_top_left_corner()};
      std::string c_tr{table_top_right_corner()};
      std::string h_7{table_horizontal_line<7>()};
      std::string h_14{table_horizontal_line<14>()};
      std::string h_23{table_horizontal_line<23>()};
      std::string h_78{table_horizontal_line<78>()};
      std::string v_ll{table_vertical_line() + " "};
      std::string v_rr{" " + table_vertical_line()};
      std::string v_lc{" " + table_vertical_line() + " "};
      std::string j_tt{table_top_junction()};
      std::string j_cc{table_center_cross()};
      std::string j_ll{table_left_junction()};
      std::string j_rr{table_right_junction()};
 
      *this->m_ostream
        << c_tl << h_78 << c_tr << std::endl
        << v_ll << "Solver:" << std::setw(69) << this->name() << v_rr << std::endl
        << v_ll << "Task:  " << std::setw(69) << this->task() << v_rr << std::endl
        << 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
        << v_ll << "#Iter" << v_lc << "   #f" << v_lc << "  #Df" << v_lc << " #DDf" << v_lc << " #Rlx" << v_lc
        << "   ║f(x)║₂  " << v_lc << "Additional notes" << std::setw(9) << v_rr << std::endl
        << 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;
    }

    void bottom()
    {
      std::string c_bl{table_bottom_left_corner()};
      std::string c_br{table_bottom_right_corner()};
      std::string h_7{table_horizontal_line<7>()};
      std::string h_14{table_horizontal_line<14>()};
      std::string h_23{table_horizontal_line<23>()};
      std::string h_78{table_horizontal_line<78>()};
      std::string v_ll{table_vertical_line() + " "};
      std::string v_rr{" " + table_vertical_line()};
      std::string j_ll{table_left_junction()};
      std::string j_rr{table_right_junction()};
      std::string j_bb{table_bottom_junction()};
 
      *this->m_ostream
        << 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
        << v_ll << std::setw(40) << (this->m_converged ? "CONVERGED" : "NOT CONVERGED") << std::setw(40) << v_rr << std::endl
        << c_bl << h_78 << c_br << std::endl;
    }

    void info(Real residuals, std::string const & notes = "-")
    {
      std::string v_rr{" " + table_vertical_line()};
      std::string v_ll{table_vertical_line() + " "};
      std::string v_lc{" " + table_vertical_line() + " "};
 
      *this->m_ostream << v_ll
        << std::setw(5) << this->m_iterations << v_lc
        << std::setw(5) << this->m_function_evaluations << v_lc
        << std::setw(5) << this->m_first_derivative_evaluations << v_lc
        << std::setw(5) << this->m_second_derivative_evaluations << v_lc
        << std::setw(5) << this->m_relaxations << v_lc
        << std::setw(12) << std::scientific << std::setprecision(6) << residuals << v_lc
        << notes << std::setw(25-notes.length()) << v_rr << std::endl;
    }
 
  }; // class Solver
 
} // namespace Optimist
 
#endif // OPTIMIST_SOLVER_HH