RootFinderDerivative.h

#ifndef theplu_yat_utility_root_finder_derivative_
#define theplu_yat_utility_root_finder_derivative_

// $Id: RootFinderDerivative.h 4243 2022-09-18 23:25:28Z peter $

/*
  Copyright (C) 2022 Peter Johansson

  This file is part of the yat library, https://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <https://www.gnu.org/licenses/>.
*/

#include "Exception.h"
#include "univariable/f.h"
#include "univariable/df.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

#include <limits>
#include <memory>
#include <sstream>

namespace theplu {
namespace yat {
namespace utility {

  class RootFinderDerivative
  {
  public:
    RootFinderDerivative(const RootFinderDerivative&) = delete;

    RootFinderDerivative& operator=(const RootFinderDerivative&) = delete;

    unsigned int epochs(void) const;

    const char* name(void) const;

    class Stopper
    {
    public:
      virtual bool operator()(const gsl_root_fdfsolver*)=0;
    };

    class Delta : public Stopper
    {
    public:
      Delta(double epsabs, double epsrel);

      bool operator()(const gsl_root_fdfsolver*);
    private:
      double epsabs_;
      double epsrel_;
      std::unique_ptr<double> prev_;
    };

    template<class FUNC>
    class Residual : public Stopper
    {
    public:
      Residual(const FUNC& func, double epsabs)
        : func_(func), epsabs_(epsabs)
      {}

      bool operator()(const gsl_root_fdfsolver*);
    private:
      const FUNC& func_;
      double epsabs_;
    };

    template<class FUNC>
    double operator()(FUNC& func, double guess, Stopper&& stopper);

    template<class FUNC>
    double operator()(FUNC& func, double guess, Stopper&& stopper,
                      unsigned int max_epochs);

  protected:
    RootFinderDerivative(const gsl_root_fdfsolver_type* t);
  private:
    const gsl_root_fdfsolver_type* type_;
    struct GslFree
    {
      void operator()(gsl_root_fdfsolver*) const;
    };
    std::unique_ptr<gsl_root_fdfsolver, GslFree> solver_;
    unsigned int epochs_;
  };


  // template implementation
  template<class FUNC>
  double
  RootFinderDerivative::operator()(FUNC& func, double guess,
                                   RootFinderDerivative::Stopper&& stopper)
  {
    unsigned int max_epochs = std::numeric_limits<unsigned int>::max();
    return (*this)(func, guess, std::move(stopper), max_epochs);
  }


  template<class FUNC>
  double
  RootFinderDerivative::operator()(FUNC& func, double guess,
                                   RootFinderDerivative::Stopper&& stopper,
                                   unsigned int max_epochs)
  {
    gsl_function_fdf gsl_func;
    gsl_func.f = univariable::f<FUNC>;
    gsl_func.df = univariable::df<FUNC>;
    gsl_func.fdf = univariable::fdf<FUNC>;
    gsl_func.params = &func;
    gsl_root_fdfsolver_set(solver_.get(), &gsl_func, guess);

    for (epochs_=0; epochs_<max_epochs && !stopper(solver_.get()); ++epochs_) {
      int status = gsl_root_fdfsolver_iterate(solver_.get());
      if (status) {
        std::ostringstream ss;
        ss << "RootFinderDerivative: at x = "
           << gsl_root_fdfsolver_root(solver_.get());
        throw GSL_error(ss.str(), status);

      }
    }

    return gsl_root_fdfsolver_root(solver_.get());
  }


  template<class FUNC>
  bool
  RootFinderDerivative::Residual<FUNC>::operator()(const gsl_root_fdfsolver* s)
  {
    double x = gsl_root_fdfsolver_root(s);
    double f = func_(x);
    return gsl_root_test_residual(f, epsabs_) == GSL_SUCCESS;
  }

}}}
#endif