RootFinder.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_ROOTFINDER_HH
#define OPTIMIST_ROOTFINDER_HH
 
#include "Optimist/SolverBase.hh"
 
namespace Optimist
{

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

    template <typename Real, Integer N, typename DerivedSolver, bool ForceEigen = false>
    class RootFinder : public SolverBase<Real, N, N, DerivedSolver, ForceEigen>
    {
    public:
      // Fancy static assertions (just for fun, don't take it too seriously)
      static_assert(N != 0, "Are you sure you want to solve a zero-dimensional system of equations?");
 
      friend class SolverBase<Real, N, N, RootFinder<Real, N, DerivedSolver, ForceEigen>>;
 
      static constexpr bool is_rootfinder{true};
      static constexpr bool is_optimizer{false};
 
      OPTIMIST_BASIC_CONSTANTS(Real)
 
      // I/O types
      using Vector = typename SolverBase<Real, N, N, DerivedSolver, ForceEigen>::InputType;
 
      // Derivative types
      using Matrix = typename SolverBase<Real, N, N, DerivedSolver, ForceEigen>::FirstDerivativeType;
      using Tensor = typename SolverBase<Real, N, N, DerivedSolver, ForceEigen>::SecondDerivativeType;

      RootFinder() {}

      std::string name() const {return static_cast<const DerivedSolver *>(this)->name_impl();}

      Integer jacobian_evaluations() const {return this->first_derivative_evaluations();}

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

      void max_jacobian_evaluations(Integer t_jacobian_evaluations)
      {
        this->max_first_derivative_evaluations(t_jacobian_evaluations);
      }

      Integer hessian_evaluations() const {return this->first_derivative_evaluations();}

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

      void max_hessian_evaluations(Integer t_hessian_evaluations)
      {
        this->max_first_derivative_evaluations(t_hessian_evaluations);
      }
 
    protected:
      template <typename JacobianLambda>
      bool evaluate_jacobian(JacobianLambda && jacobian, Vector const & x, Matrix & out)
      {
        return this->evaluate_first_derivative(std::forward<JacobianLambda>(jacobian), x, out);
      }

      template <typename HessianLambda>
      bool evaluate_hessian(HessianLambda && hessian, Vector const & x, Matrix & out)
      {
        return this->evaluate_second_derivative(std::forward<HessianLambda>(hessian), x, out);
      }
 
    public:
      template <typename FunctionLambda>
      bool solve(FunctionLambda && function, Vector const & x_ini, Vector & x_sol)
      {
        #define CMD "Optimist::RootFinder::solve(...): "
 
        static_assert(DerivedSolver::requires_function,
          CMD "the solver requires a function.");
        return static_cast<DerivedSolver *>(this)->solve_impl(
          std::forward<FunctionLambda>(function),
          x_ini, x_sol);
 
        #undef CMD
      }

      template <typename FunctionLambda, typename JacobianLambda>
      bool solve(FunctionLambda && function, JacobianLambda && jacobian, Vector const & x_ini,
        Vector & x_sol)
      {
        #define CMD "Optimist::RootFinder::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_impl(
          std::forward<FunctionLambda>(function),
          std::forward<JacobianLambda>(jacobian),
          x_ini, x_sol);
 
        #undef CMD
      }

      template <typename FunctionLambda, typename JacobianLambda, typename HessianLambda>
      bool solve(FunctionLambda && function, JacobianLambda && jacobian, HessianLambda && hessian,
        Vector const & x_ini, Vector & x_sol)
      {
        #define CMD "Optimist::RootFinder::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_impl(
          std::forward<FunctionLambda>(function),
          std::forward<JacobianLambda>(jacobian),
          std::forward<HessianLambda>(hessian),
          x_ini, x_sol);
 
        #undef CMD
      }
 
    }; // class RootFinder
 
    // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    template <typename Real, typename DerivedSolver>
    class RootFinder<Real, 1, DerivedSolver> : public SolverBase<Real, 1, 1, DerivedSolver>
    {
    public:
      friend class SolverBase<Real, 1, 1, RootFinder<Real, 1, DerivedSolver>>;
 
      static constexpr bool is_rootfinder{true};
      static constexpr bool is_optimizer{false};
 
      OPTIMIST_BASIC_CONSTANTS(Real)
 
      
      RootFinder<Real, 1, DerivedSolver>() {}

      std::string name() const {return static_cast<const DerivedSolver *>(this)->name_impl();}

      template <typename FunctionLambda>
      bool solve(FunctionLambda && function, Real x_ini, Real & x_sol)
      {
        return static_cast<DerivedSolver *>(this)->solve_impl(
          std::forward<FunctionLambda>(function),
          x_ini, x_sol);
      }

       template <typename FunctionLambda, typename FirstDerivativeLambda>
      bool solve(FunctionLambda && function, FirstDerivativeLambda && first_derivative, Real x_ini,
        Real & x_sol)
      {
        return static_cast<DerivedSolver *>(this)->solve_impl(
          std::forward<FunctionLambda>(function),
          std::forward<FirstDerivativeLambda>(first_derivative),
          x_ini, x_sol);
      }

      template <typename FunctionLambda, typename FirstDerivativeLambda, typename SecondDerivativeLambda>
      bool solve(FunctionLambda && function, FirstDerivativeLambda && first_derivative, SecondDerivativeLambda
        && second_derivate, Real x_ini, Real & x_sol)
      {
        return static_cast<DerivedSolver *>(this)->solve_impl(
          std::forward<FunctionLambda>(function),
          std::forward<FirstDerivativeLambda>(first_derivative),
          std::forward<SecondDerivativeLambda>(second_derivate),
          x_ini, x_sol);
      }
 
    }; // class RootFinder
 
  } // namespace RootFinder
 
} // namespace Optimist
 
#endif // OPTIMIST_ROOTFINDER_HH