Utils_AlgoBracket.hh

Go to the documentation of this file.

/*--------------------------------------------------------------------------*\
 |                                                                          |
 |  Copyright (C) 2025                                                      |
 |                                                                          |
 |         , __                 , __                                        |
 |        /|/  \               /|/  \                                       |
 |         | __/ _   ,_         | __/ _   ,_                                |
 |         |   \|/  /  |  |   | |   \|/  /  |  |   |                        |
 |         |(__/|__/   |_/ \_/|/|(__/|__/   |_/ \_/|/                       |
 |                           /|                   /|                        |
 |                           \|                   \|                        |
 |                                                                          |
 |      Enrico Bertolazzi                                                   |
 |      Dipartimento di Ingegneria Industriale                              |
 |      Università degli Studi di Trento                                    |
 |      email: enrico.bertolazzi@unitn.it                                   |
 |                                                                          |
\*--------------------------------------------------------------------------*/
 
//
// file: Utils_AlgoBracket.hh
//
 
#pragma once
 
#ifndef UTILS_ALGO_BRACKET_dot_HH
#define UTILS_ALGO_BRACKET_dot_HH
 
#include <iostream>
#include <sstream>
#include <string>
#include <algorithm>
#include <cmath>
 
#include "Utils.hh"
 
namespace Utils {
 
  using std::pow;
  using std::abs;
  using std::sqrt;

 
 
  /*
  //   ____                 _        _
  //  | __ ) _ __ __ _  ___| | _____| |_
  //  |  _ \| '__/ _` |/ __| |/ / _ \ __|
  //  | |_) | | | (_| | (__|   <  __/ |_
  //  |____/|_|  \__,_|\___|_|\_\___|\__|
  //
  */

  template <typename Real>
  class Bracket_base_fun {
  public:
    virtual Real eval( Real x ) const = 0;
    Real operator () ( Real x ) const { return this->eval(x); }
  };
 
  #ifndef DOXYGEN_SHOULD_SKIP_THIS
  template <typename Real, typename PFUN>
  class Bracket_fun : public Bracket_base_fun<Real> {
    PFUN m_fun;
  public:
    explicit Bracket_fun( PFUN pfun ) : m_fun(pfun) {}
    Real eval( Real x ) const override { return m_fun(x); };
  };
  #endif

  template <typename Real>
  class AlgoBracket {
  
    using Method = enum class Method : unsigned {
      BISECTION=0, ILLINOIS, CHANDRUPATLA, BRENT, RIDDER, MODIFIED_AB, ALGO748
    };
 
    using Integer = int;
 
    Real m_tolerance_x{ 100*machine_eps<Real>() };
    Real m_tolerance_f{ 100*machine_eps<Real>() };
    bool m_converged{ false };
 
    Real   m_a{0}, m_fa{0};
    Real   m_b{0}, m_fb{0};
    Method m_select{Method::CHANDRUPATLA};
 
    Bracket_base_fun<Real> * m_function{nullptr};
 
    Integer m_max_fun_evaluation{1000}; // max number of function evaluations
    Integer m_max_iteration{200};       // max number of iterations
 
    mutable Integer m_iteration_count{0};    // explore iteration counter
    mutable Integer m_fun_evaluation_count{0};
 
    Real evaluate( Real x ) { ++m_fun_evaluation_count; return m_function->eval(x); };
 
    Real eval();
    Real eval( Real a, Real b );
    Real eval( Real a, Real b, Real amin, Real bmax );
 
    Real p_zero2( Real d, Real fd ) const;
    Real invp_zero2( Real d, Real fd ) const;
    Real invp_zero3( Real d, Real fd, Real e, Real fe ) const;
 
  public:
 
    AlgoBracket() = default;
    ~AlgoBracket() = default;
 
    explicit
    AlgoBracket( Real tol_x, Real tol_f )
    : m_tolerance_x( tol_x ), m_tolerance_f( tol_f )
    {}

    Real
    eval( Real a, Real b, Bracket_base_fun<Real> * fun ) {
      m_function = fun;
      return this->eval( a, b );
    }

    Real
    eval( Real a, Real b, Real amin, Real bmax, Bracket_base_fun<Real> * fun ) {
      m_function = fun;
      return this->eval( a, b, amin, bmax );
    }

    template <typename PFUN>
    Real
    eval2( Real a, Real b, PFUN pfun ) {
      Bracket_fun<Real,PFUN> fun( pfun );
      m_function = &fun;
      return this->eval( a, b );
    }

    template <typename PFUN>
    Real
    eval2( Real a, Real b, Real amin, Real bmax, PFUN pfun ) {
      Bracket_fun<Real,PFUN> fun( pfun );
      m_function = &fun;
      return this->eval( a, b, amin, bmax );
    }

    template <typename PFUN>
    Real
    eval3( Real a, Real b, Real fa, Real fb, PFUN pfun ) {
      Bracket_fun<Real,PFUN> fun( pfun );
      m_function             = &fun;
      m_iteration_count      = 0;
      m_fun_evaluation_count = 0;
      m_a = a; m_fa = fa;
      m_b = b; m_fb = fb;
      return eval();
    }

    void set_max_iterations( Integer mit );

    void set_max_fun_evaluation( Integer mfev );

    Integer used_iter() const { return m_iteration_count; }

    Integer num_fun_eval() const { return m_fun_evaluation_count; }

    Real tolerance_x() const { return m_tolerance_x; }

    Real tolerance_f() const { return m_tolerance_f; }

    bool converged() const { return m_converged; }
 
    Real a()  const { return m_a; }
    Real b()  const { return m_b; }
    Real fa() const { return m_fa; }
    Real fb() const { return m_fb; }
 
    string algo() const {
      switch ( m_select ) {
      case Method::BISECTION:    return "bisection";    break;
      case Method::ILLINOIS:     return "illinois";     break;
      case Method::CHANDRUPATLA: return "Chandrupatla"; break;
      case Method::BRENT:        return "Brent";        break;
      case Method::RIDDER:       return "Ridder";       break;
      case Method::MODIFIED_AB:  return "modified_AB";  break;
      case Method::ALGO748:      return "algo748";      break;
      }
      return "";
    }
 
    void select( unsigned i_algo ) {
      switch ( i_algo ) {
      case 0: m_select = Method::BISECTION;    break;
      case 1: m_select = Method::ILLINOIS;     break;
      case 2: m_select = Method::CHANDRUPATLA; break;
      case 3: m_select = Method::BRENT;        break;
      case 4: m_select = Method::RIDDER;       break;
      case 5: m_select = Method::MODIFIED_AB;  break;
      case 6: m_select = Method::ALGO748;      break;
      }
    }
 
  };
 
  #ifndef UTILS_OS_WINDOWS
  extern template class AlgoBracket<float>;
  extern template class AlgoBracket<double>;
  #endif

 
}
 
#endif
 
//
// EOF: Utils_AlgoBracket.hh
//