linesearch.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
#include <kernel/base_header.hpp>
#include <kernel/solver/iterative.hpp>
 
#include <deque>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Functional_, typename Filter_>
    class Linesearch : public IterativeSolver<typename Functional_::VectorTypeR>
    {
      public:
        typedef Filter_ FilterType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename Functional_::DataType DataType;
        typedef IterativeSolver<typename Functional_::VectorTypeR> BaseClass;
 
      protected:
        // Note that this cannot be const, as the functional saves its state and thus changes
        Functional_& _functional;
        Filter_& _filter;
 
        VectorType _vec_grad;
        VectorType _vec_initial_sol;
        VectorType _vec_tmp;
        VectorType _vec_pn;
 
        DataType _fval_min;
        DataType _fval_0;
        DataType _trim_threshold;
 
        DataType _alpha_0;
        DataType _alpha_min;
        DataType _delta_0;
        DataType _norm_dir;
        DataType _norm_sol;
 
        DataType _tol_curvature;
        DataType _tol_decrease;
        DataType _tol_step;
        bool _dir_scaling;
 
      public:
        std::deque<VectorType>* iterates;
 
      public:
        explicit Linesearch(const String& plot_name, Functional_& functional, Filter_& filter,
          bool keep_iterates = false) :
          BaseClass(plot_name),
          _functional(functional),
          _filter(filter),
          _fval_min(Math::huge<DataType>()),
          _fval_0(Math::huge<DataType>()),
          _trim_threshold(Math::huge<DataType>()),
          _alpha_0(1),
          _alpha_min(0),
          _delta_0(Math::huge<DataType>()),
          _norm_dir(0),
          _norm_sol(0),
          _tol_curvature(DataType(0.3)),
          _tol_decrease(DataType(1e-3)),
          _tol_step(DataType(Math::pow(Math::eps<DataType>(), DataType(0.85)))),
          _dir_scaling(false),
          iterates(nullptr)
          {
            // set communicator by functional (same interface as matrix)
            this->_set_comm_by_matrix(functional);
 
            if(keep_iterates)
            {
              iterates = new std::deque<VectorType>;
            }
          }
 
        explicit Linesearch(const String& plot_name, const String& section_name, const PropertyMap* section,
        Functional_& functional, Filter_& filter) :
          BaseClass(plot_name, section_name, section),
          _functional(functional),
          _filter(filter),
          _fval_min(Math::huge<DataType>()),
          _fval_0(Math::huge<DataType>()),
          _trim_threshold(Math::huge<DataType>()),
          _alpha_0(1),
          _alpha_min(DataType(0)),
          _delta_0(Math::huge<DataType>()),
          _norm_dir(DataType(0)),
          _norm_sol(DataType(0)),
          _tol_curvature(0.3),
          _tol_decrease(1e-3),
          _tol_step(DataType(Math::pow(Math::eps<DataType>(), DataType(0.85)))),
          _dir_scaling(false),
          iterates(nullptr)
          {
            // Check if we have to keep the iterates
            auto keep_iterates_p = section->query("keep_iterates");
            if(keep_iterates_p.second)
            {
              iterates = new std::deque<VectorType>;
            }
 
            auto tol_curvature_p = section->query("tol_curvature");
            if(tol_curvature_p.second)
            {
              set_tol_curvature(DataType(std::stod(tol_curvature_p.first)));
            }
 
            auto tol_decrease_p = section->query("tol_decrease");
            if(tol_decrease_p.second)
            {
              set_tol_decrease(DataType(std::stod(tol_decrease_p.first)));
            }
            // Check if we have to keep the iterates
            auto tol_step_p = section->query("tol_step");
            if(tol_step_p.second)
            {
              set_tol_step(DataType(std::stod(tol_step_p.first)));
            }
          }
 
        virtual ~Linesearch()
        {
          if(iterates != nullptr)
          {
            delete iterates;
          }
        }
 
        virtual String get_summary() const override
        {
          String msg(this->get_plot_name()+ ": its: "+stringify(this->get_num_iter())
              +" ("+ stringify(this->get_status())+")"
              +", evals: "+stringify(_functional.get_num_func_evals())+" (func) "
              + stringify(_functional.get_num_grad_evals()) + " (grad) "
              + stringify(_functional.get_num_hess_evals()) + " (hess)"
              +" last step: "+stringify_fp_sci(_alpha_min)+"\n");
          msg +=this->get_plot_name()+": fval: "+stringify_fp_sci(_fval_0)
            + " -> "+stringify_fp_sci(_fval_min)
            + ", factor "+stringify_fp_sci(_fval_min/_fval_0)+"\n";
          msg += this->get_plot_name()  +": <dir, grad>: "+stringify_fp_sci(this->_def_init)
            + " -> "+stringify_fp_sci(this->_def_cur)
            + ", factor " +stringify_fp_sci(this->_def_cur/this->_def_init);
          return msg;
        }
 
        void set_tol_curvature(DataType tol_curvature)
        {
          XASSERT(tol_curvature > DataType(0));
 
          _tol_curvature = tol_curvature;
        }
 
        void set_tol_decrease(DataType tol_decrease)
        {
          XASSERT(tol_decrease > DataType(0));
 
          _tol_decrease = tol_decrease;
        }
 
        void set_tol_step(DataType tol_step)
        {
          XASSERT(tol_step > DataType(0));
 
          _tol_step = tol_step;
        }
 
        void set_dir_scaling(const bool b)
        {
          _dir_scaling = b;
        }
 
        virtual void init_symbolic() override
        {
          // Create temporary vectors
          _vec_initial_sol = this->_functional.create_vector_r();
          _vec_tmp = this->_functional.create_vector_r();
          _vec_pn = this->_functional.create_vector_r();
          _vec_grad = this->_functional.create_vector_r();
          BaseClass::init_symbolic();
        }
 
        virtual void done_symbolic() override
        {
          // Clear temporary vectors
          _vec_initial_sol.clear();
          _vec_pn.clear();
          _vec_tmp.clear();
          _vec_grad.clear();
          BaseClass::done_symbolic();
        }
 
        virtual void reset()
        {
          _fval_min = Math::huge<DataType>();
          _fval_0 = Math::huge<DataType>();
          _trim_threshold = Math::huge<DataType>();
          _alpha_0 = DataType(1),
          _alpha_min = DataType(0);
          _delta_0 = Math::huge<DataType>();
          _norm_dir = DataType(0);
          _norm_sol = DataType(0);
 
          if(iterates != nullptr)
          {
            iterates->clear();
          }
        }
 
        virtual DataType get_rel_update() const
        {
          return Math::abs(this->_alpha_min);
        }
 
        DataType get_final_fval() const
        {
          return _fval_min;
        }
 
        void set_initial_fval(DataType f0)
        {
          _fval_0 = f0;
          if(_trim_threshold == Math::huge<DataType>())
          {
            _trim_threshold = DataType(10)*(Math::abs(_fval_0) + DataType(1));
          }
        }
 
        void get_defect_from_grad(VectorType& vec_def) const
        {
          vec_def.scale(this->_vec_grad, DataType(-1));
        }
 
        void set_grad_from_defect(const VectorType& vec_def)
        {
          this->_vec_grad.scale(vec_def, DataType(-1));
        }
 
        const VectorType& get_initial_sol() const
        {
          return _vec_initial_sol;
        }
 
        virtual void trim_func_grad(DataType& func)
        {
          if(func > _trim_threshold)
          {
            func = _trim_threshold;
            this->_vec_grad.format(DataType(0));
          }
        }
 
      protected:
 
        virtual Status _startup(DataType& alpha, DataType& fval, DataType& delta, const VectorType& vec_sol, const VectorType& vec_dir)
        {
          Status status(Status::progress);
 
          this->_num_iter = Index(0);
          this->_vec_initial_sol.copy(vec_sol);
 
          this->_vec_pn.copy(vec_dir);
          this->_norm_dir = this->_vec_pn.norm2();
 
          // First scale so that all entries are |.| < 1
          this->_vec_pn.scale(this->_vec_pn, DataType(1)/this->_vec_pn.max_abs_element());
          // Now scale so that vec_pn.norm2() == 1. Note that this will only be approximately true
          this->_vec_pn.scale(this->_vec_pn, DataType(1)/this->_vec_pn.norm2());
 
          // It is critical that fval_0 was set from outside using set_initial_fval!
          this->_fval_min = this->_fval_0;
          this->_delta_0 = this->_vec_pn.dot(this->_vec_grad);
 
          if(this->_delta_0 > DataType(0))
          {
            XABORTM("Search direction is not a descent direction: " +stringify_fp_sci(this->_delta_0));
          }
 
          // Compute initial defect. We want to minimize d^T * grad(_functional)
          this->_def_init = Math::abs(this->_delta_0);
 
          // Norm of the initial guess
          this->_norm_sol = vec_sol.norm2();
 
          this->_alpha_min = DataType(0);
 
          // plot?
          if(this->_plot_iter())
          {
            String msg = this->_plot_name
              +  ": " + stringify(this->_num_iter).pad_front(this->_iter_digits)
              + " : " + stringify_fp_sci(this->_def_init)
              + " : " + stringify_fp_sci(this->_fval_0)
              + ", ||dir|| = " + stringify_fp_sci(this->_norm_dir);
            this->_print_line(msg);
          }
 
          if(!Math::isfinite(this->_fval_0) || !Math::isfinite(this->_delta_0) || !Math::isfinite(this->_norm_dir))
          {
            status = Status::aborted;
          }
 
          // Set initial step size
          alpha = _alpha_0;
          // Other intitial values
          fval = _fval_0;
          delta = _delta_0;
 
          return status;
        }
 
        virtual Status _check_convergence(const DataType fval, const DataType df, const DataType alpha)
        {
          Status status(Status::progress);
 
          this->_def_cur = Math::abs(df);
 
          Statistics::add_solver_expression(
            std::make_shared<ExpressionDefect>(this->name(), this->_def_cur, this->get_num_iter()));
 
          // plot?
          if(this->_plot_iter())
          {
            String msg = this->_plot_name
              +  ": " + stringify(this->_num_iter).pad_front(this->_iter_digits)
              + " : " + stringify_fp_sci(this->_def_cur)
              + " / " + stringify_fp_sci(this->_def_cur / this->_def_init)
              + " : " + stringify_fp_sci(fval) + " : " + stringify_fp_sci(alpha);
            this->_print_line(msg);
          }
 
          // ensure that the defect is neither NaN nor infinity
          if(!Math::isfinite(this->_def_cur))
          {
            status = Status::aborted;
          }
 
          // is diverged?
          if(this->is_diverged())
          {
            status = Status::diverged;
          }
 
          // If the maximum number of iterations was performed, return the iterate for the best step so far
          if(this->_num_iter >= this->_max_iter)
          {
            status = Status::max_iter;
          }
 
          // If the strong Wolfe conditions hold, we are successful
          if((fval < this->_fval_0 +_tol_decrease*alpha*_delta_0) && (Math::abs(df) < -_tol_curvature*_delta_0))
          {
            status = Status::success;
          }
 
          return status;
        }
 
    }; // class Linesearch
  } // namespace Solver
} // namespace FEAT